-- test.lua -- Typed 2016-01-23 by H. Henkel -- Public domain mpfr = require("lmpfrlib") mpfr.set_default_prec(10000) mpfr.set_default_rounding_mode("MPFR_RNDA") local a = mpfr.num(2) local b = mpfr.num(3) print("================================") print(a) print("================================") if a == b then print("=== equal") else print("=== not equal") end c = a + b mpfr.printf("%.50Rf\n", c) d = -c mpfr.printf("%.50Rf\n", d) e = c - d mpfr.printf("%.50Rf\n", e) f = e * a mpfr.printf("%.50Rf\n", f) g = f / b mpfr.printf("%.50Rf\n", g) h = g ^ g mpfr.printf("%.50Rf\n", h) i = h ^ 3 mpfr.printf("%.50Rf\n", i) j = mpfr.const_euler() print("euler:") mpfr.printf("%.50Rf\n", j) k = mpfr.const_log2() print("log2:") mpfr.printf("%.50Rf\n", k) ------------------------------------------------------------------------ local EK = function(k) -- A&S 17.6 local pi = mpfr.const_pi() local k_prime = mpfr.sqrt(1 - k^2) -- A&S 17.2.18, 17.2.19 local a, b, c = {}, {}, {} a[0], b[0], c[0] = mpfr.num(1), k_prime, k -- A&S 17.6.2 local n local csum = c[0]^2 -- A&S 17.6.4 for i = 1, 100 do n = i a[i] = 0.5 * (a[i - 1] + b[i - 1]) -- A&S 17.6.1 b[i] = mpfr.sqrt(a[i - 1] * b[i - 1]) c[i] = 0.5 * (a[i - 1] - b[i - 1]) csum = csum + c[i]^2 * 2^i if c[i] < mpfr.num("1e-20000") then break end end assert(n < 100) local K = pi / (2 * a[n]) -- A&S 17.6.3 local E = K * (1 - 0.5 * csum) -- A&S 17.6.4 return E, K end E, K = EK(mpfr.num("0.999999999999999999999999999999999999999999")) print("E=" .. mpfr.snprintf(1000, "%.70Rf\n", E)) print("K=" .. mpfr.snprintf(1000, "%.70Rf\n", K))