Matroids Matheplanet Forum Index
Moderiert von matroid
Mathematik » Numerik & Optimierung » Fast and robust calculation of Jacobi theta functions with complex arguments
Autor
Kein bestimmter Bereich Fast and robust calculation of Jacobi theta functions with complex arguments
MontyPythagoras
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 13.05.2014
Mitteilungen: 2901
Wohnort: Werne
  Themenstart: 2022-05-02 12:43

Dear all, I am writing this post in English because I have contacted Fredrik Johansson, a French mathematician who is the original author of the mpmath library for the Python programming language, facilitating arbitrary-precision floating-point arithmetic, and I would like him to be able understand at least this initial post. Everyone else please feel free to contribute in German. I found in Fredriks github repository the unanswered question if there were ways to efficiently calculate the Jacobi theta functions with the nome $q$ approaching 1. Recently we had a thread here on the matheplanet peripherically touching that matter. Since there is a variety of differing definitions, I am hereafter using the definition given on Wolfram mathworld, equations (2) through (5), or (6) through (11) respectively, while the latter set of equations matches the definitions in the Digital library of mathematical functions. When $|q|\to1$ and $\Im(z)\neq0$, the Jacobi theta function quickly exceeds the limitations of double precision floats. Therefore, the calculation of the natural logarithm of the theta function allows to easily surpass the double precision limit. The general goal is to apply transformations to the arguments of the function in order to move $q$ towards smaller absolute values, thus accelerating convergence. Also, it is beneficial to use a slightly different notation than in the references, reducing the need for multiplications and divisions by $\pi$. Hence, with $q=e^{i\pi\tau}$ and $z=\pi\zeta$ I will be using the following definition: $$\tau\equiv\frac{\ln q}{i\pi}\tag1$$$$\zeta\equiv\frac z\pi\tag2$$$$\vartheta_3(z,q)\equiv\vartheta_3(\zeta,\tau)=\sum_{n=-\infty}^{\infty}e^{i\pi(n^2\tau+2n\zeta)}\tag3$$The larger $\Im(\tau)$, the faster the convergence. Also, it is important to know that even though there are four definitions of seemingly different functions, each of them can be expressed through any of the others. Since $\vartheta_3$ is the easiest to code, the first step is to express $\vartheta_{1,2,4}$ as transformations of $\vartheta_3$, as follows: $$\vartheta_1(\zeta,\tau)=\exp\left(i\pi\left(\zeta+\frac\tau4-\frac12\right)\right)\cdot\vartheta_3\left(\zeta+\frac\tau2-\frac12,\tau\right)\tag4$$$$\vartheta_2(\zeta,\tau)=\exp\left(i\pi\left(\zeta+\frac\tau4\right)\right)\cdot\vartheta_3\left(\zeta+\frac\tau2,\tau\right)\tag5$$$$\vartheta_4(\zeta,\tau)=\vartheta_3\left(\zeta-\frac12,\tau\right)\tag6$$(see here and consider $\vartheta_3(\zeta,\tau)=\vartheta_3(\zeta-1,\tau)$). A transformation of arguments can generally be expressed as $$\vartheta_3(\zeta_k,\tau_k)=f_k(\zeta_k,\tau_k)\cdot \vartheta_3(\zeta_{k+1},\tau_{k+1})\tag7$$with $f_k$ just being a pre-factor that is dependent on arguments $\zeta_k$ and $\tau_k$, with $k$ representing the $k^{\text{th}}$ transformation step. Apart from the initial transformation in equations (4) through (6), there are three important transformations to perform repeatedly until a certain criterion is fulfilled. The first one is $$\vartheta_3(\zeta,\tau)=\exp\left(i\pi\left(n^2\tau-2n\zeta\right)\right)\cdot\vartheta_3(\zeta-m-n\tau,\tau)\tag8$$derived from equation 20.2.8 here. This transformation allows to keep $\Im(\zeta)$ small-ish, and it is easily performed in two steps: a) $n=\Big\lfloor\frac{\Im(\zeta)}{\Im(\tau)}\Big\rfloor$ b) $m=\lfloor\Re(\zeta)\rfloor$ Since $\tau$ will likely have a real part as well, it is important to perform b) after a) is completed. The second important transformation is $$\vartheta_3(\zeta,\tau)=\sqrt{-i\tau'}\exp\left(i\pi\zeta^2\tau'\right)\cdot\vartheta_3(\zeta\tau',\tau')=\sqrt{\tau'}\exp\left(i\pi\left(\zeta^2\tau'-\frac14\right)\right)\cdot\vartheta_3(\zeta\tau',\tau')\tag9$$with $$\tau'=-\frac1\tau\tag{10}$$(equation 20.7.32 here). At first glance it may seem pointless to perform this transformation twice in a row, since apparently $\tau''=\tau$, but that is not the case, because after performing this transformation once on a $\tau$ that has both a small real and a small imaginary part, $\tau'$ will have a larger imaginary part, but the larger real part of $\tau'$ can be diminished by subtracting or adding any even natural number, always keeping the real portion small, and after the next transformation the imaginary part is even larger, and so forth. Lastly, to further improve transformation (9), we use $$\vartheta_3(\zeta,\tau)=\vartheta_3\left(\zeta-\frac12,\tau\pm1\right)\tag{11}$$The following picture illustrates, how consecutive application of transformations (10) and (11) move the nome $q$ towards zero: The original nome is $q_0=0.59999+0.8i$, after seven transformations it is $q_7=0.0010036351+0.0052840721i$ with $|q|=0.0053785408$. Various tests have shown that convergence happens sufficiently fast when $\Im(\tau)>0.4$ (equivalent to $|q|<0.2846...$). As I mentioned above, I am calculating the natural logarithm of the theta function. Hence, (7) can be written as $$\ln\vartheta_3(\zeta_0,\tau_0)=\sum_{k=0}^{n-1}\ln f_k(\zeta_k,\tau_k)+ \ln\vartheta_3(\zeta_{n},\tau_{n})\tag{12}$$Most of the pre-factors above had the form $\exp(i\pi\cdot(\dots))$, the logarithm of which can easily be calculated without using time-consuming $\exp$ or $\log$ functions. Only the first factor ($\sqrt{\tau'}$) in equation (9) would need to be logarithmized. It is easier and faster to use two separate pre-factors, so we rewrite (7) as $$\vartheta_3(\zeta_k,\tau_k)=f_k(\zeta_k,\tau_k)\cdot\exp\left(i\pi g_k(\zeta_k,\tau_k)\right)\cdot \vartheta_3(\zeta_{k+1},\tau_{k+1})\tag{13}$$Mathematically, this is futile, but the code gets a lot faster by avoiding $\exp$ or $\log$ whenever possible. With this re-definition we can write generally $$\ln\vartheta_3(\zeta_0,\tau_0)=i\pi\sum_{k=0}^{n-1}g_k(\zeta_k,\tau_k)+ \ln\left(\vartheta_3(\zeta_{n},\tau_{n})\prod_{k=0}^{n-1}f_k(\zeta_k,\tau_k)\right)\tag{14}$$Finally, after having reached a favorable set of arguments, we still need to calculate $\vartheta_3(\zeta_n,\tau_n)$ by using (3). We can directly determine the lower and upper limit of the summation indices. With $\tau=\tau_r+i\tau_i$ and $\zeta=\zeta_r+i\zeta_i$, (3) can be written as $$\vartheta_3(\zeta,\tau)=\sum_{n=n_l}^{n_u}\exp\left(i\pi(n^2\left(\tau_r+i\tau_i)+2n(\zeta_r+i\zeta_i)\right)\right)=\sum_{n=n_{low}}^{n_{up}}\exp\left(-\pi(n^2\tau_i+2n\zeta_i)\right)\cdot\exp\left(i\pi(n^2\tau_r+2n\zeta_r)\right)\tag{15}$$The second exponential function in this sum has an absolute value of 1 and is irrelevant for convergence. The first exponential function on the other hand determines the absolute value of the summand. If the summand is to be of any significance for the entire sum, we can demand that $$\exp\left(-\pi(n^2\tau_i+2n\zeta_i)\right)>2^{-c}\tag{16}$$For double precision we can set $c=55$ because the mantissa of a double precision float only has 53 bits, and if (16) holds true, any summands outside of the range between $(n_{low}\dots n_{up})$ can be neglected. We will show that $n=0$ will always be within this range, so that there is always a summand with an absolute value of 1. From (16) we get $$-(n^2\tau_i+2n\zeta_i)>-\frac{c\ln2}{\pi}$$$$n_{low}=\Big\lceil-\frac{\zeta_i}{\tau_i}-\sqrt{\left(\frac{\zeta_i}{\tau_i}\right)^2+\frac{c\ln2}{\pi\tau_i}}\Big\rceil=-\Big\lfloor\frac{\zeta_i}{\tau_i}+\sqrt{\left(\frac{\zeta_i}{\tau_i}\right)^2+\frac{c\ln2}{\pi\tau_i}}\Big\rfloor\leq0\tag{17}$$$$n_{up}=\Big\lfloor-\frac{\zeta_i}{\tau_i}+\sqrt{\left(\frac{\zeta_i}{\tau_i}\right)^2+\frac{c\ln2}{\pi\tau_i}}\;\Big\rfloor\geq0\tag{18}$$(Keep in mind that $\tau_i=\Im(\tau)>0$ and therefore, the square root is larger than the absolute value of $\frac{\zeta_i}{\tau_i}$). With this set of equations, I programmed the Jacobi theta functions in Python and I would like to share and dissect (and maybe improve) the code here. Okay, here goes: \sourceon Python \numberson import math, cmath PI = math.pi PIj = complex(0, math.pi) TAU_1 = 1 / math.tau TAUj = complex(0, math.tau) def ln_theta(n, z, q): if q == 0: if n < 3: return -math.inf else: return 0 if abs(q) >= 1: return math.nan tau = cmath.log(q) / PIj z /= PI f = 1 if n == 1 or n == 4: z -= 0.5 if n < 3: g = z + 0.25 * tau z += 0.5 * tau else: g = 0 n = z.imag // tau.imag if n != 0: a = n * tau g += n * (a - 2 * z) z -= a z -= int(z.real) while tau.imag < 0.4: n = round(tau.real) if n & 1: z -= 0.5 tau = -1 / (tau - n) g += tau * z * z - 0.25 f *= cmath.sqrt(tau) z *= tau q, n = math.modf(z.imag / tau.imag) if n != 0: a = n * tau g += n * (a - 2 * z) z -= a z -= int(z.real) s = 0 z *= 2 a = math.sqrt(q * q + 12.135 / tau.imag) # 12.135 ≈ 55 * ln(2) / PI for n in range(-int(a + q), int(a - q) + 1): s += cmath.exp(PIj * n * (tau * n + z)) s = PIj * g + cmath.log(f * s) return s - TAUj * round(TAU_1 * s.imag) \sourceoff 1. Lines 9 to 14: filtering out trivial cases 2. Lines 15 and 16: equations (1) and (2) 3. Lines 17 to 23: initialization and equations (4) through (6) 4. Lines 25 to 29: equation 8 with $n$ according to a), line 26 only saving time 5. Line 30: equation 8 with $m$ according to b) 6. Line 32: "while" loop until $\Im(\tau)>0.4$ (values above 0.35 don't seem to make a difference, values below 0.35 increase the duration) 7. Lines 33 to 35: equations (10) and (11) 8. Lines 36 to 38: equation (9) 9. Lines 40 to 45: same as 25 to 30, but at the same time storing the remainder of $\frac{\Im(\zeta)}{\Im(\tau)}$ for later calculation of the summation limits in variable $q$ 10. Lines 47 to 50: equation (15), limits using (17) and (18). The value 12,135 in line 49 applies to double precision floats only and would need to be increased in case of higher precision targets 11. Line 51: equation (14) 12. Line 52: normalization of the imaginary part of the logarithm to $(-\pi,\pi)$ and return On my computer with an 8th generation i7 processor, the calculation takes under 7µs for "harmless" arguments and below 12µs for critical ones, in a regular Python interpreter. In a pre-compiled code like C++ or CPython, calculation might very well finish in nano seconds. Except for line 15 at the beginning and lines 50 and 51 at the end, complex $\exp$ and $\log$ functions could be avoided, the algorithm basically only needs elementary arithmetic and square roots. I would be interested in comparisons regarding speed and accuracy. Ciao, Thomas


   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1677
  Beitrag No.1, eingetragen 2022-05-02 20:00

Hallo Thomas, interessant wird es ja dann, wenn man die Nullstellen von Gleichungen 5. Grades per QuinticEquation und den Thetha-Funktionen sofort ausrechnen kann. Jedoch haben die Gleichungen (7)...(12) eine Hilfsvariable p und q, die ich momentan nur in der Sonderfallgleichung (13) sehe. Oder gibt es Transformationen oder bessere Quellen, um wirklich universell alle Faktoren a0...a5 einsetzen zu können? Ansonsten gibt es ja für Spezialfälle die einfacheren hypergeometrischen Funktionen: https://www.lamprechts.de/gerd/Bilder/PolynomLoesungen/xH5_16.png Grüße Gerd


   Profil
MontyPythagoras
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 13.05.2014
Mitteilungen: 2901
Wohnort: Werne
  Beitrag No.2, vom Themenstarter, eingetragen 2022-05-14 11:14

Hallo Gerd, \quoteon(2022-05-02 20:00 - hyperG in Beitrag No. 1) Hallo Thomas, interessant wird es ja dann, wenn man die Nullstellen von Gleichungen 5. Grades per QuinticEquation und den Thetha-Funktionen sofort ausrechnen kann. \quoteoff die Theta-Funktionen haben vielfältige Anwendungen, nicht nur die quintischen Gleichungen. Darum soll es hier aber nicht gehen. Ich habe den Code angepasst, so dass ich mithilfe der mpmath library die Theta-Funktionen mit beliebiger Präzision relativ schnell berechnen kann. Hier ist der Code (aus Platzspargründen ausgeblendet): \showon \sourceon Python \numberson from mpmath import sqrt, pi, log, exp, expjpi, nan, mpc, fabs, frac, mp, nint, jtheta, nstr from time import process_time as prtm def theta(n, z, q): if q == 0: return 0 if n < 3 else 1 if fabs(q) >= 1: return nan currprec = mp.prec crit = log(2) / pi * (currprec + 2) mp.prec += 8 + round(2.3 * log((1 + z.imag * z.imag) / (1 - fabs(q)))) tau = -1j * log(q) / pi z /= pi f = 1 if n == 1 or n == 4: z -= 0.5 if n < 3: g = z + 0.25 * tau z += 0.5 * tau else: g = 0 n = nint(z.imag / tau.imag) if n != 0: a = z z -= n * tau g -= n * (a + z) while tau.imag < 0.5: n = nint(tau.real) z -= frac(n * 0.5) z -= nint(z.real) tau = -1 / (tau - n) a = z z *= tau g += a * z - 0.25 f *= sqrt(tau) q = z.imag / tau.imag n = nint(q) if n != 0: q -= n a = z z -= n * tau g -= n * (a + z) s = 0 z *= 2 a = sqrt(q * q + crit / tau.imag) for n in range(-int(a + q), int(a - q) + 1): s += expjpi(n * (tau * n + z)) s *= f * expjpi(g) mp.prec = currprec return s # End of function definition def looptime(f, a, t_test): l = 0 t = prtm() while prtm() - t < 0.5: r = f(* a) l += 1 t = prtm() - t if t < 0.5 * t_test: l = int(t_test / t * l) + 1 t = prtm() for k in range(l): r = f(* a) t = prtm() - t return r, t / l def shortnum(x, left, right): a = nstr(x, mp.dps, strip_zeros=False) if a.find("j") > -1: l1 = a.find(" ") a, b = a[:l1], a[l1:] l1 = a.find(".") + 1 + left l2 = (max(0, a.find("e"))-right) % len(a) r = a if l2 - l1 < 3 else a[:l1] + "_(" + str(l2 - l1) + ")_" + a[l2:] l1 = b.find(".") + 1 + left l2 = (max(0, b.find("e"))-right) % (len(b) - 2) r = r + b if l2 - l1 < 3 else r + b[:l1] + "_(" + str(l2 - l1) + ")_" + b[l2:] else: l1 = a.find(".") + 1 + left l2 = (max(0, a.find("e"))-right) % len(a) r = a if l2 - l1 < 3 else a[:l1] + "_(" + str(l2 - l1) + ")_" + a[l2:] return r mp.dps = 200 t_test = 5 #seconds test duration target a = 2, mpc("0.5", "1"), mpc("0.599999", "0.8") print("n = {:d}\nz = {}\nq = {}".format(a[0], nstr(a[1], 15), nstr(a[2], 15))) print("Precision: {} digits\n".format(mp.dps)) r, t_2 = looptime(theta, a, t_test) print(" theta = {} in {:6.2f} ms".format(shortnum(r, 4, 8), t_2 * 1000)) r, t_1 = looptime(jtheta, a, t_test) print("jtheta = {} in {:6.2f} ms".format(shortnum(r, 4, 8), t_1 * 1000)) print("\nFactor {:.2f}".format(t_1 / t_2)) \sourceoff \showoff Die Berechnung für diese Funktionsargumente mit $|q|\to1$ ist mit einer Berechnungsdauer von unter 5 Millisekunden über 3000 mal schneller als die in mpmath implementierte Funktion "jtheta". Hier ist der Output des obigen Codes: \sourceon output n = 2 z = (0.5 + 1.0j) q = (0.599999 + 0.8j) Precision: 200 digits theta = (2.9177_(187)_30027968e+723825 - 6.4828_(187)_82484441e+723825j) in 4.89 ms jtheta = (2.9177_(187)_30027968e+723825 - 6.4828_(187)_82484441e+723825j) in 14812.50 ms Factor 3027.48 \sourceoff Rückt man mit $q$ noch näher an 1, erzeugt "jtheta" einen Fehler. Auch Wolframalpha streckt die Glieder, wenn man sich mehr Stellen anzeigen lassen will. Ich bin gespannt, was Mathematica schafft. Machen wir es also noch fieser (die implementierte jtheta versagt hier den Dienst): \sourceon output n = 2 z = (0.5 + 1.0j) q = (0.59999999 + 0.8j) Precision: 1000 digits theta = (2.3477_(987)_65680478e+72382415 + 1.7577_(987)_64645963e+72382415j) in 78.87 ms \sourceoff Knapp 79ms für 1000 Stellen Genauigkeit bei $q = 0.59999999 + 0.8i$. Nicht schlecht für eine Interpreter-Sprache und ein komplexes Ergebnis mit über 72 Millionen Stellen vor dem Komma... 🙃 (Die Funktion shortnum kürzt die Ausgabe ab. "_(987)_" bedeutet, dass 987 Stellen ausgeblendet wurden.) Ciao, Thomas


   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1677
  Beitrag No.3, eingetragen 2022-05-14 15:33

Da hast Du ja einen "schönen Argumenten-Spezialfall" gefunden, der richtig schlecht konvergiert! Mathematica 12 ist bei 200 Stellen nur bis 5999/10^4 zu gebrauchen: \sourceon mathematica Block[{$MaxExtraPrecision=\[Infinity]},N[EllipticTheta[2,1/2+1*I,5999/10^4+8/10*I],200]]//Timing Out[1]= 4.140625 s 2.2004985806624864352...54188597044*10^7239 +2.714583554422487342...89360687584*10^7239 I \sourceoff Schon 59999/10^5 geht in den Minuten-Bereich! Was mach eigentlich Deine Funktion expjpi(x) ? Ah, hab's gefunden bei https://mpmath.org/doc/current/functions/powers.html: Exp[x*I*Pi]


   Profil
MontyPythagoras
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 13.05.2014
Mitteilungen: 2901
Wohnort: Werne
  Beitrag No.4, vom Themenstarter, eingetragen 2022-05-14 17:33

Hallo Gerd, \quoteon(2022-05-14 15:33 - hyperG in Beitrag No. 3) Da hast Du ja einen "schönen Argumenten-Spezialfall" gefunden \quoteoff Nö, den muss man nicht finden. Das klingt nach Zufall, aber die Konvergenz wird grundsätzlich immer sehr schleppend, wenn sich der Betrag von $q$ dem Wert 1 nähert. Ich wollte halt einen komplexen Wert für $q$, der knapp kleiner als 1 ist, und da bot sich das pythagoreische Tripel $(0.6, 0.8, 1)$ an. Jeder andere Wert knapp unter 1 ist also genauso schlimm, und ein imaginärer Anteil in $z$ macht das Ganze noch viel schlimmer. Mit den von Dir benutzten Werten erhalte ich: \sourceon output n = 2 z = (0.5 + 1.0j) q = (0.5999 + 0.8j) Precision: 200 digits theta = (2.200498580662486435_(170)_54185386122e+7239 + 2.714583554422487342_(170)_89362373610e+7239j) in 5.04 ms jtheta = (2.200498580662486435_(170)_54185386122e+7239 + 2.714583554422487342_(170)_89362373610e+7239j) in 1515.62 ms Factor 300.54 \sourceoff Hier ist mein Code "nur" 300 mal schneller als jtheta der mpmath Bibliothek. Der Vorteil wird um so größer, je näher $|q|$ an 1 rückt. Aber immerhin über 800 mal schneller als Mathematica. Was allerdings interessant ist, ist die Tatsache, dass mein Algo und die jtheta Funktion von mpmath exakt die gleichen Werte ausspucken, die auch stimmen, wenn man eine deutlich größere Genauigkeit eingibt und dann die Stellen um die Position 200 herum vergleicht. Mathematica wird auf den letzten 7 Stellen ungenau. Würde mich interessieren, wie sich die anderen Programme wie Maple etc. schlagen. Ciao, Thomas


   Profil
MontyPythagoras
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 13.05.2014
Mitteilungen: 2901
Wohnort: Werne
  Beitrag No.5, vom Themenstarter, eingetragen 2022-05-14 17:56

Korrektur: bei nochmaligen Vergleich stimmen alle Endziffern von Mathematica, die von Python weichen ab. Gebe ich allerdings eine "precision" von 207 ein, stimmt das Ergebnis mit Mathematica in den ersten 200 Stellen überein, bei einer Berechnungsdauer von 5,10ms. Ciao, Thomas


   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1677
  Beitrag No.6, eingetragen 2022-05-14 18:06

Frage zu int(a + q) Wenn a & q komplex sind, was soll denn bei int für ein Ganzzahlwert rauskommen? [Die Antwort wurde nach Beitrag No.3 begonnen.]


   Profil
MontyPythagoras
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 13.05.2014
Mitteilungen: 2901
Wohnort: Werne
  Beitrag No.7, vom Themenstarter, eingetragen 2022-05-14 20:49

Hallo Gerd, $a$ und $q$ sind reelle Zahlen, siehe oben Gleichungen (17) und (18). Ciao, Thomas


   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1677
  Beitrag No.8, eingetragen 2022-05-14 22:21

\quoteon(2022-05-14 20:49 - MontyPythagoras in Beitrag No. 7) Hallo Gerd, $a$ und $q$ sind reelle Zahlen, siehe oben Gleichungen (17) und (18). Ciao, Thomas \quoteoff Hi Thomas, klar -> es war ein ganz anderer Fehler bei mir: statt "If" hatte ich an einer stelle "IF" und schon zog sich der Folgefehler durch alle Zeilen und brachte alles durcheinander :-) https://matheplanet.com/matheplanet/nuke/html/images/forum/subject/rotate.gif Nun funktioniert es (letzte Stellen hatte ich noch nicht genau analysiert, da ich zunächst mp.prec += 8 + round(2.3 * log((1 + z.imag * z.imag) / (1 - fabs(q)))) weggelassen hatte) \showon https://matheplanet.com/matheplanet/nuke/html/uploads/b/47407_Mtheta.png \showoff Nun "kann es Mathematica auch schneller" für diese speziellen Argumente :-) Interessanter Algorithmus. Für n=3 hatte ich ihn noch nicht getestet. Grüße Gerd


   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1677
  Beitrag No.9, eingetragen 2022-05-14 23:13

Nun zur Genauigkeit: Um mit anderen Argumenten (QuinticEquation mit 2. Argument=0) auf 50 Stellen genau zu sein, muss man intern mit 131 Stellen rechnen: \sourceon mathematica Mtheta[3, 0, 59995/10^5 + 7/10*I, 131] - Block[{$MaxExtraPrecision = [Infinity]}, N[EllipticTheta[3, 0, 59995/10^5 + 7/10*I], 50]] 10^-50 \sourceoff Das ist natürlich gewaltig (mehr als das Doppelte!) und mit 3 Argumenten recht kompliziert vorherzubestimmen...


   Profil
MontyPythagoras
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 13.05.2014
Mitteilungen: 2901
Wohnort: Werne
  Beitrag No.10, vom Themenstarter, eingetragen 2022-05-15 16:38

Hallo Gerd, da scheint noch ein Fehler in Deiner Umsetzung des Algorithmus zu sein. Vergleiche aus Deinem Post #8 das Ergebnis für $q=0.5999+0.8i$ mit dem aus Deinem Beitrag #3. Den Wert aus Beitrag #3 kann ich bestätigen. Auch die anderen beiden aus #8 stimmen nicht komplett. Das könnte der Grund sein, dass Du intern mit so vielen Stellen rechnen musst. Das kann ich nämlich auch nicht bestätigen. Ich schaue mir bei Gelegenheit Deinen Algorithmus noch einmal an, bin allerdings kein Mathematica-Experte. Ciao, Thomas


   Profil
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1677
  Beitrag No.11, eingetragen 2022-05-15 21:52

\quoteon(2022-05-15 16:38 - MontyPythagoras in Beitrag No. 10) ...anderen beiden aus #8 stimmen nicht komplett. ... \quoteoff Das hatte ich doch beschrieben, dass ich "MontyGenau": mp.prec +=8+round(2.3*log((1+z.imag *z.imag)/(1-fabs(q)))) absichtlich weggelassen hatte, um die wirklich nötige interne Genauigkeit exakt bestimmen zu können. Du hattest bisher immer nur mit n=2 gerechnet! Mit steigendem n muss intern noch weiter angehoben werden, was diese Tabelle zeigt: \sourceon nameDerSprache n = 2; z = 1/2 + 1*I; q = 5999/10^4 + 8/10*I; genau = 50; noetigGenau = 150; MontyGenau=82 n = 3; z = 0; q = 59995/10^5 + 7/10*I; genau = 50; noetigGenau = 131; MontyGenau=64 n = 4; z = 0; q = 59995/10^5 + 79/100*I;genau = 50; noetigGenau = 155; MontyGenau=69 \sourceoff Was natürlich sein kann, ist die Tatsache, dass Dein crit = log(2) / pi * (currprec + 2) was ich ja 1:1 übernommen hatte, nicht mit meiner unveränderten Genauigkeit zusammen passt, da Du ja den Code für MontyGenau optimiert hattest. Ich vermute, dass es so lauten müsste, damit wir die identische Anzeige bekommen wollen, wenn ich intern 1:1 rechne {und den Genauigkeitsparameter "außen" vergrößert übergebe}: crit = log(2) / pi * (MontyGenau + 2) Fakt ist aber, dass n bei der Berechnung der wirklich benötigten internen Genauigkeit mit berücksichtigt werden muss. Wie gesagt: bei 3 Argumenten ist das nicht einfach. Grüße Gerd


   Profil
MontyPythagoras
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 13.05.2014
Mitteilungen: 2901
Wohnort: Werne
  Beitrag No.12, vom Themenstarter, eingetragen 2022-05-15 22:42

Hallo Gerd, ich glaube, da hast Du einiges falsch verstanden. Das $n$, welches als Argument übergeben wird, identifiziert den Typ der Theta-Funktion, denn es gibt derer 4, siehe Gleichungen (3) bis (6). Die erste Transformation besteht darin, $\vartheta_1$, $\vartheta_2$ und $\vartheta_4$ durch $\vartheta_3$ auszudrücken, weil letzteres am einfachsten zu programmieren ist. Wir berechnen also nur $\vartheta_3$, nach einer passenden Transformation der Argumente! Ob Du also hier mit $n=2$ oder $n=1, 3, 4$ testest, ist völlig egal. Nach dieser ersten Transformation ist der Variablenbuchstabe $n$ freigeworden. Ich hätte, um Verwechslungen dieser Art zu vermeiden, einen anderen Variablenbuchstaben wählen können, aber ich bin faul und eine weitere Variable hätte auch zusätzlichen Speicher belegt. Daher bekommt $n$ nach Zeile 19 in Beitrag #2 eine komplett neue Bedeutung. Didaktisch unklug, ich weiß, aber das $n=2$ aus den Startargumenten ist ab Zeile 19 komplett unbedeutend und hat mit den Genauigkeiten überhaupt nichts zu tun. Erst in Zeile 48 in #2 wird eigentlich die übriggebliebene Funktion $\vartheta_3$ berechnet, und dort werden durch Gleichungen (17) und (18) die Grenzen der Reihe berechnet, also welche Summanden mit einbezogen werden müssen . DAS ist entscheidend für die Genauigkeit. Dort liegt bei Dir schon mal ein kleiner Fehler vor, nämlich bei dem "IntegerPart[a-q]+1". Das "+1" am Ende ist unnötig, es resultiert im Python-Code aus der Tatsache, dass ein "range" immer eine Nummer vorher endet. In Python umfasst der "range(2,5)" die Elemente 2, 3 und 4, aber die 5 nicht. Ich glaube aber, ich weiß, wo das Problem liegt. In der Python-Library mpmath gibt es zwei Möglichkeiten, die Genauigkeit zu setzen, und zwar "mp.dps" und "mp.prec". Der Unterschied ist, dass "dps" die Anzahl der Dezimalstellen betrifft (hier also 200), "prec" dagegen die Anzahl der Bits. Ich nutze intern in dem Algorithmus also Bits statt Dezimalstellen. Du müsstest bei Dir, wenn Du in Dezimalstellen denkst, crit=log(10)/Pi*... statt crit=log(2)/Pi*... setzen. 10 statt 2 im Logarithmus. Dann wird es passen. Ciao, Thomas


   Profil
MontyPythagoras hat die Antworten auf ihre/seine Frage gesehen.

Wechsel in ein anderes Forum:
 Suchen    
 
All logos and trademarks in this site are property of their respective owner. The comments are property of their posters, all the rest © 2001-2022 by Matroids Matheplanet
This web site was originally made with PHP-Nuke, a former web portal system written in PHP that seems no longer to be maintained nor supported. PHP-Nuke is Free Software released under the GNU/GPL license.
Ich distanziere mich von rechtswidrigen oder anstößigen Inhalten, die sich trotz aufmerksamer Prüfung hinter hier verwendeten Links verbergen mögen.
Lesen Sie die Nutzungsbedingungen, die Distanzierung, die Datenschutzerklärung und das Impressum.
[Seitenanfang]