-- b2-mpfr.lua -- -- Formula for the Inductance -- of a Helix with Wire of any Section -- Chester Snow, Physicist, Bureau of Standards -- November 10, 1926 -- Washington, Government Printing Office, 1926 -- -- Here: Calculation of constant B2 (90) on p. 458 -- Test case for Lua MPFR API, allocation-efficient version -- -- Written 2016-01-23, changed 2016-01-27 -- Copyright (C) 2016 Hartmut Henkel -- E-Mail: hartmut_henkel@gmx.de -- Homepage: http://www.hhenkel.de/lmpfrlib/ -- -- This program is free software: you can redistribute it and/or modify -- it under the terms of the GNU Lesser Public License as published by -- the Free Software Foundation, either version 3 of the License, or -- (at your option) any later version. -- -- This program is distributed in the hope that it will be useful, -- but WITHOUT ANY WARRANTY; without even the implied warranty of -- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -- GNU Lesser Public License for more details. -- -- You should have received a copy of the GNU Lesser Public License -- along with this program. If not, see . mpfr = require("lmpfrlib") mpfr.set_default_prec(1000) mpfr.set_default_rounding_mode("MPFR_RNDA") local B2 = function() local zero = mpfr.init_set_si(0) local onehalf = mpfr.init_set_d(0.5) local one = mpfr.init_set_si(1) -- 4 * ln(2) local fourln2 = mpfr.init_set_si(2) fourln2:log(fourln2) fourln2:mul_si(fourln2, 4) -- local B2 = 0 local B2 = mpfr.init_set(zero) -- n = 1 local n = mpfr.init_set(one) -- local T0_nm1 = 4 * ln(2) local T0_nm1 = mpfr.init_set(fourln2) -- local sum_nm1 = 1 / (n - 1 + 1 / 2) local sum_nm1 = mpfr.init() local denom = mpfr.init() denom:sub_si(n, 1) denom:add(denom, onehalf) sum_nm1:si_div(1, denom) -- local T0_n = 4 * ln(2) - 2 * sum_nm1 local T0_n = mpfr.init_set(fourln2) T0_n:sub(T0_n, sum_nm1) T0_n:sub(T0_n, sum_nm1) -- local sum_n = sum_nm1 + 1 / (n + 1 / 2) local sum_n = mpfr.init() denom:add(n, onehalf) local frac = mpfr.init() frac:si_div(1, denom) sum_n:add(sum_nm1, frac) -- local T0_np1 = 4 * ln(2) - 2 * sum_n local T0_np1 = mpfr.init_set(fourln2) T0_np1:sub(T0_np1, sum_n) T0_np1:sub(T0_np1, sum_n) -- local B2 = B2 + (T0_nm1 + T0_np1) / n^2 local num = mpfr.init() num:add(T0_nm1, T0_np1) local denom = mpfr.init() denom:sqr(n) frac:div(num, denom) B2:add(B2, frac) local j = mpfr.init_set(zero) local res = mpfr.init() local loop = mpfr.init_set_si(1000000) while(true) do n:add_si(n, 1) j:add_si(j, 1) -- T0_nm1 = T0_n T0_nm1:set(T0_n) -- T0_n = T0_np1 T0_n:set(T0_np1) -- sum_nm1 = sum_n sum_nm1:set(sum_n) -- sum_n = sum_nm1 + 1 / (n + 1 / 2) denom:add(n, onehalf) frac:si_div(1, denom) sum_n:add(sum_nm1, frac) -- T0_np1 = 4 * ln(2) - 2 * sum_n T0_np1:set(fourln2) T0_np1:sub(T0_np1, sum_n) T0_np1:sub(T0_np1, sum_n) -- B2 = B2 + (T0_nm1 + T0_np1) / n^2 num:add(T0_nm1, T0_np1) denom:sqr(n) frac:div(num, denom) B2:add(B2, frac) -- if j == loop then print("# " .. mpfr.snprintf(1000, "%-50.0Rf", n)) res:div_d(B2, 8) local s = mpfr.snprintf(1000, "%.50Rf", res) print(s) j:set(zero) end end end B2()