-- lether-wenston-mpfr.lua -- -- Abscissas and Weights for Gauss-Legendre Quadrature -- -- Zeros calculated by (L&W): -- Minimax approximations to the zeros of Pn(x) -- and Gauss-Legendre quadrature -- F. G. Lether, P. R. Wenston -- Journal of Computational and Applied Mathematics 59 (1995) 245--252 -- -- Weights calculated by (D&R): -- Abscissas and Weights for Gaussian Quadratures of High Order -- P. Davis and P. Rabinowitz -- Journal of Research of the National Bureau of Standards -- Vol. 56, No. 1, January 1956, Research Paper 2645 -- -- Written 2016-01-23, Version 2017-06-29 -- Copyright (C) 2016, 2017 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 . ------------------------------------------------------------------------ local version = "2017-06-29" local m = require("lmpfrlib") local sqrt, cos, pi = m.sqrt, m.cos, m.const_pi() local one = m.num(1) -- L&W p. 245 local P = function(x, n) assert(type(x) == "userdata") assert(type(n) == "number") local Pmm local Pm = one -- P0(x) local P = x -- P1(x) if n == 0 then return Pm elseif n == 1 then return P else for i = 2, n do Pmm, Pm = Pm, P P = (1 - one / i) * (x * Pm - Pmm) + x * Pm end end return P, Pm end -- L&W p. 245 local P_faster = function(x, n) assert(type(x) == "userdata") assert(type(n) == "number") local Pmm = m.init() local Pm = m.init() local P = m.init() local aux1 = m.init() local aux2 = m.init() local aux3 = m.init() Pm:set(one) -- P0(x) P:set(x) -- P1(x) if n == 0 then return Pm elseif n == 1 then return P else for i = 2, n do Pmm:set(Pm) Pm:set(P) -- P = (1 - one / i) * (x * Pm - Pmm) + x * Pm aux1:div_si(one, i) aux1:neg(aux1) aux1:add_si(aux1, 1) aux2:mul(x, Pm) aux3:set(aux2) aux2:sub(aux2, Pmm) aux1:mul(aux1, aux2) P:add(aux1, aux3) end end return P, Pm end -- L&W p. 247 local E = function(x, n) assert(type(x) == "userdata") assert(type(n) == "number") local P, Pm = P_faster(x, n) local v = P / (n * (Pm - x * P)) local B = one / 3 * (3 + n * (n + one)) * x^2 - one / 3 * (1 + n * (n + one)) local C = one / 6 * (6 + 5 * n * (n + one)) * x^3 - one / 6 * (4 + 5 * n * (n + one)) * x assert(type(v) == "userdata") assert(type(B) == "userdata") assert(type(C) == "userdata") local E = x - (1 - x) * (1 + x) * v * (1 + v * (x + v * (B + C * v))) return E, Pm end -- L&W p. 246 local Theta = function(n, k) assert(type(n) == "number") assert(type(k) == "number") assert((k >= 1) and (k <= n // 2)) local Theta = (k - one / 4) / (n + one / 2) * pi -- (1.4) return Theta end -- L&W p. 251 local xn1 = function(n) assert(type(n) == "number") local j1 = m.num("2.4048255577") local a1 = m.num("0.0105077313672") local b1 = m.num("0.165291040565") local xn1 = cos(j1 / sqrt(n^2 + n + one / 3) * (1 - a1 / (360 * (n^2 + n + b1)^2))) -- Example 5. return xn1 end -- L&W p. 247 local Omega = function(n, k) assert(type(n) == "number") assert(type(k) == "number") local Theta = Theta(n, k) local n = m.num(n) local Omega = (1 - one / 8 * n^(-2) + m.num(5) / 38 * n^(-3) - m.num(2) / 25 * n^(-4) * (1 - m.num(14) / 39 * Theta^(-2))) * cos(Theta) -- (1.5) return Omega end local Gauss = function(n, k) local x, Pm if k == 1 then x = xn1(n) else x = Omega(n, k) end for i = 1, 3 do x, Pm = E(x, n) -- (1.6) end local w = 2 * (1 - x^2) / (n * Pm)^2 -- D&R (4) return x, w end ------------------------------------------------------------------------ local gen_lua_dqg = function(n, prec, mpfrprec) m.set_default_prec(mpfrprec) m.set_default_rounding_mode("MPFR_RNDA") assert((n // 2) == (n / 2)) local x, w = {}, {} for i = 1, n // 2 do x[i], w[i] = Gauss(n, i) end local f = io.open("dqg" .. n .. ".lua", "w") f:write("-- dqg" .. n .. ".lua\n") f:write("-- Abscissas and Weights for Gauss-Legendre Quadrature\n") f:write("-- This file is generated by lether-wenston-mpfr.lua\n") f:write("-- Version " .. version .. "\n") f:write("-- MPFR default precision = " .. mpfrprec .. " bit\n") f:write("\n") f:write("-- abscissas\n") f:write("local x = {\n") for i = 1, n // 2 do f:write(" " .. m.sprintf("%." .. prec .. "Rf", x[i]) .. ",\n") end f:write("}\n") f:write("\n") f:write("-- weights\n") f:write("local w = {\n") for i = 1, n // 2 do f:write(" " .. m.sprintf("%." .. prec .. "Rf", w[i]) .. ",\n") end f:write("}\n") f:write("\n") f:write("local dqg" .. n .. " = function (xl, xu, fct)\n") f:write(" local a = 0.5 * (xu + xl)\n") f:write(" local b = 0.5 * (xu - xl)\n") f:write(" local y = 0\n") f:write(" for i = 1, " .. n // 2 .. " do\n") f:write(" local c = b * x[i]\n") f:write(" y = y + w[i] * (fct(a + c) + fct(a - c))\n") f:write(" end\n") f:write(" y = y * b\n") f:write(" return y\n") f:write("end\n") f:write("\n") f:write("return dqg" .. n .. "\n") f:close() end ------------------------------------------------------------------------ local gen_fortran_dqg = function(n, prec, mpfrprec) m.set_default_prec(mpfrprec) m.set_default_rounding_mode("MPFR_RNDA") assert((n // 2) == (n / 2)) local x, w = {}, {} for i = 1, n // 2 do x[i], w[i] = Gauss(n, i) end local f = io.open("dqg" .. n .. ".f", "w") local subname = "DQG" .. n local subcall = "DQG" .. n .. " (XL,XU,FCT,Y)" local subrout = "DQG" .. n .. "(XL,XU,FCT,Y)" local i = 10 f:write(string.format("C DG%06d\n", i)) i = i + 10 f:write(string.format("C ..................................................................DG%06d\n", i)) i = i + 10 f:write(string.format("C DG%06d\n", i)) i = i + 10 f:write(string.format("C SUBROUTINE %-25s DG%06d\n", subname, i)) i = i + 10 f:write(string.format("C DG%06d\n", i)) i = i + 10 f:write(string.format("C PURPOSE DG%06d\n", i)) i = i + 10 f:write(string.format("C TO COMPUTE INTEGRAL(FCT(X), SUMMED OVER X FROM XL TO XU) DG%06d\n", i)) i = i + 10 f:write(string.format("C DG%06d\n", i)) i = i + 10 f:write(string.format("C USAGE DG%06d\n", i)) i = i + 10 f:write(string.format("C CALL %-25s DG%06d\n", subcall, i)) i = i + 10 f:write(string.format("C PARAMETER FCT REQUIRES AN EXTERNAL STATEMENT DG%06d\n", i)) i = i + 10 f:write(string.format("C DG%06d\n", i)) i = i + 10 f:write(string.format("C DESCRIPTION OF PARAMETERS DG%06d\n", i)) i = i + 10 f:write(string.format("C XL - DOUBLE PRECISION LOWER BOUND OF THE INTERVAL. DG%06d\n", i)) i = i + 10 f:write(string.format("C XU - DOUBLE PRECISION UPPER BOUND OF THE INTERVAL. DG%06d\n", i)) i = i + 10 f:write(string.format("C FCT - THE NAME OF AN EXTERNAL DOUBLE PRECISION FUNCTION DG%06d\n", i)) i = i + 10 f:write(string.format("C SUBPROGRAM USED. DG%06d\n", i)) i = i + 10 f:write(string.format("C Y - THE RESULTING DOUBLE PRECISION INTEGRAL VALUE. DG%06d\n", i)) i = i + 10 f:write(string.format("C DG%06d\n", i)) i = i + 10 f:write(string.format("C REMARKS DG%06d\n", i)) i = i + 10 f:write(string.format("C NONE DG%06d\n", i)) i = i + 10 f:write(string.format("C DG%06d\n", i)) i = i + 10 f:write(string.format("C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED DG%06d\n", i)) i = i + 10 f:write(string.format("C THE EXTERNAL DOUBLE PRECISION FUNCTION SUBPROGRAM FCT(X) DG%06d\n", i)) i = i + 10 f:write(string.format("C MUST BE FURNISHED BY THE USER. DG%06d\n", i)) i = i + 10 f:write(string.format("C DG%06d\n", i)) i = i + 10 f:write(string.format("C METHOD DG%06d\n", i)) i = i + 10 local str = n .. "-POINT GAUSS QUADRATURE" f:write(string.format("C EVALUATION IS DONE BY MEANS OF %-29sDG%06d\n", str, i)) i = i + 10 local str = "" .. 2 * n - 1 f:write(string.format("C FORMULA, WHICH INTEGRATES POLYNOMIALS UP TO DEGREE %-9sDG%06d\n", str, i)) i = i + 10 f:write(string.format("C EXACTLY. FOR REFERENCE, SEE DG%06d\n", i)) i = i + 10 f:write(string.format("C V.I.KRYLOV, APPROXIMATE CALCULATION OF INTEGRALS, DG%06d\n", i)) i = i + 10 f:write(string.format("C MACMILLAN, NEW YORK/LONDON, 1962, PP.100-111 AND 337-340. DG%06d\n", i)) i = i + 10 f:write(string.format("C DG%06d\n", i)) i = i + 10 f:write(string.format("C COLOPHON DG%06d\n", i)) i = i + 10 f:write(string.format("C THIS FILE IS GENERATED BY LETHER-WENSTON-MPFR.LUA DG%06d\n", i)) i = i + 10 f:write(string.format("C VERSION %-52sDG%06d\n", version, i)) i = i + 10 local str = "" .. mpfrprec .. " BIT" f:write(string.format("C MPFR DEFAULT PRECISION = %-35sDG%06d\n", str, i)) i = i + 10 f:write(string.format("C ..................................................................DG%06d\n", i)) i = i + 10 f:write(string.format("C DG%06d\n", i)) i = i + 10 f:write(string.format(" SUBROUTINE %-25s DG%06d\n", subrout, i)) i = i + 10 f:write(string.format("C DG%06d\n", i)) i = i + 10 f:write(string.format("C DG%06d\n", i)) i = i + 10 f:write(string.format(" DOUBLE PRECISION XL,XU,Y,A,B,C,FCT DG%06d\n", i)) i = i + 10 f:write(string.format("C DG%06d\n", i)) i = i + 10 f:write(string.format(" A=.5D0*(XU+XL) DG%06d\n", i)) i = i + 10 f:write(string.format(" B=XU-XL DG%06d\n", i)) i = i + 10 local str = m.sprintf("%." .. prec .. "Rg", x[1] / 2) .. "*B" str = str.gsub(str, "e", "D") f:write(string.format(" C=%-59s DG%06d\n", str, i)) i = i + 10 local str = m.sprintf("%." .. prec .. "Rg", w[1] / 2) .. "*(FCT(A+C)+FCT(A-C))" str = str.gsub(str, "e", "D") f:write(string.format(" Y=%-59s DG%06d\n", str, i)) i = i + 10 for j = 2, (n // 2) - 1 do local str = m.sprintf("%." .. prec .. "Rg", x[j] / 2) .. "*B" str = str.gsub(str, "e", "D") f:write(string.format(" C=%-59s DG%06d\n", str, i)) i = i + 10 local str = m.sprintf("%." .. prec .. "Rg", w[j] / 2) .. "*(FCT(A+C)+FCT(A-C))" str = str.gsub(str, "e", "D") f:write(string.format(" Y=Y+%-59s DG%06d\n", str, i)) i = i + 10 end local str = m.sprintf("%." .. prec .. "Rg", x[n // 2] / 2) .. "*B" str = str.gsub(str, "e", "D") f:write(string.format(" C=%-59s DG%06d\n", str, i)) i = i + 10 local str = m.sprintf("%." .. prec .. "Rg", w[n // 2] / 2) .. "*(FCT(A+C)+FCT(A-C)))" str = str.gsub(str, "e", "D") f:write(string.format(" Y=B*(Y+%-59sDG%06d\n", str, i)) i = i + 10 f:write(string.format(" RETURN DG%06d\n", i)) i = i + 10 f:write(string.format(" END DG%06d\n", i)) i = i + 10 f:close() end local t = { E = E, xn1 = xn1, Omega = Omega, P = P_faster, gauss = gauss, gen_lua_dqg = gen_lua_dqg, gen_fortran_dqg = gen_fortran_dqg, version = version, } return t