Calculate Gaussian points and weights at arbitrary precision with the Fourier-Newton method in MPFR enabled GNU Awk

glatwgt.awk

# gawk -M -v PREC=100 -v lgaqd=1 -v nlat=120 -f glatwgt.awk

function abs(x)
{
  return (x > 0) ? x : -x
}

function gamma(y,     c)
{
  c[1] =          1 /           12
  c[2] =          1 /          288
  c[3] =       -139 /        51840
  c[4] =       -571 /      2488320
  c[5] =     163879 /    209018880
  c[6] =    5246819 /  75246796800
  c[7] = -534703531 / 902961561600
#  return 1 + (c[1] + (c[2] + (c[3] + c[4] * y) * y) * y) * y
  return 1 + (c[1] + (c[2] + (c[3] + (c[4] + (c[5] + (c[6] + c[7] * y) * y) * y) * y) * y) * y) * y
}

function ann_gamma(n,     pi, y, yh, a, g)
{
  pi = atan2(0, -1) 
  y = 1 / n
  yh = 0.5 / n
  gr = 1 / gamma(y)
  a = sqrt(2 * (2 * n + 1) / (n * pi))
  return a * gamma(yh) * (gr * gr)
}

function calc_ank(n, ank,     nh, i, k, k2, n2, k1)
# computes the Fourier coefficients of Pn as described in Swarztrauber 2002
{
  nh = n / 2
  ank[nh] = sqrt(2.0)
  for (i = 2; i <= n; i++) {
# Swarztrauber 2002
    ank[nh] *= sqrt(1 - 1/(4 * i * i)) 
# Modified
#    n2 = i * 2
#    ank[nh] *= sqrt((n2 - 1) * (n2 + 1)) / n2 
#    ank[nh] *= sqrt(((n2 - 1) / n2) * ((n2 + 1) / n2)) 
  }
  if (lgamma) ank[nh] = ann_gamma(n)
  for (k = 1; k <=nh; k++) {
#    k2 = 2 * k
    k1 = 2 * k - 1
# three below are equivalent
#    ank[nh - k] = (k2 - 1) * (2 * n - k2 + 2) / (k2 * (2 * n - k2 + 1)) * ank[nh - k + 1]
    ank[nh - k] = k1 * (n - k + 1) / (k * (2 * n - k1)) * ank[nh - k + 1]
#    ank[nh - k] = ((k1 * (n - k + 1)) / (k * (2 * n - k1))) * ank[nh - k + 1]
# below causes larger error
#    ank[nh - k] = (k1 / k) * ((n - k + 1) / (2 * n - k1)) * ank[nh - k + 1]
  }
  if (n == nh * 2) {
    ank[0] = 0.5 * ank[0]
  }
}

function cpdp(n, cp, dcp,     t1, t2, t3, t4, ncp, i, j)
# computes the Fourier coefficients of the Legendre polynomial
# p_n^0 and its derivative.
# n: degree. must be even
# cp: coefficients. 0..ncp = (n+1)/2
# dcp: coefficients. 1..ncp
{
  t1 = -1
  t2 = n + 1
  t3 = 0
  t4 = n + n + 1
  ncp = int((n+1)/2)
  cp[ncp] = 1
  for (i = 0; i < ncp; i++) {
    j = ncp - i
    t1++; t1++
    t2--
    t3++
    t4--; t4--
    cp[j-1] = (t1*t2)/(t3*t4) * cp[j]
  }
  for (i = 1; i <= ncp; i++) {
    dcp[i] = (i + i) * cp[i]
  }
}

function tpdp(n, theta, cp, dcp, pb, dpb,      cdt, sdt, kdo, cth, sth, chh, i) 
{
#  computes pn(theta) and its derivative dpb(theta)

  cdt = cos(theta + theta)
  sdt = sin(theta + theta)

  kdo = int(n/2)
  pb[0] = 0.5 * cp[0]
  dpb[0] = 0
  cth = cdt
  sth = sdt
  for (i = 1; i <= kdo; i++) {
#   pb += cp[k] * cos(2 * k * theta)
    pb[0] += cp[i] * cth
#   dpb += -2 * i * cp[k] * sin(2 * k * theta)
    dpb[0] += -dcp[i] * sth 
    chh = cdt * cth - sdt * sth
    sth = sdt * cth + cdt * sth
    cth = chh
  }
}

function dzeps(x,     a, b, c, eps)
# estimate unit roundoff in quatities of size x
{
  a = 4 / 3
  b = a - 1
  c = b + b + b
  eps = abs(c - 1)
  return eps * abs(x)
}

function gaqd(nlat, theta, wts,     eps, pi, pis2, nhalf, zero, i, nix, zlast, zhold, zprev, ddpb, sum)
{
# computes Gaussian points and weights using Fourier-Newton method
# following Swarztrauber 2002.

  eps = dzeps(1)
#  print "machine delta = ", eps
  eps = sqrt(eps)
  eps = eps * sqrt(eps)

  pi = atan2(0, -1) 
  pis2 = 0.5 * pi
  nhalf = int(nlat/2)
#  print "pi =", pi, " pis2 =", pis2, " nhalf = ", nhalf

  cpdp(nlat, cp, dcp)

# estimate first point next to pi/2

  zero = pis2 - 0.5 * (pi / nlat)
#  print "zero =", zero
  for (i = 0; i < nhalf; i++) {
    nix = nhalf - i - 1

# Newton iteration

    zlast = zero
    tpdp(nlat, zero, cp, dcp, pb, dpb)
    zero += -pb[0]/dpb[0]
    while (abs(zero - zlast) > eps * abs(zero)) {
      zlast = zero
      tpdp(nlat, zero, cp, dcp, pb, dpb)
      zero += -pb[0]/dpb[0]
    }
    theta[nix] = zero
    zhold = zero

# Yakimiw's formula permits using old pb and dpb

    ddpb = dpb[0] + pb[0] * cos(zlast)/sin(zlast)
    wts[nix] = (nlat + nlat + 1) / (ddpb * ddpb)
    if (i == 0) {
      zero = 3 * zero - pi
    } else {
      zero = zero + zero - zprev
    }
   zprev = zhold
 }

# Extend points and weights via symmetries

  for (i = 0; i < nhalf; i++) {
    wts[nlat - i] = wts[i]
    theta[nlat -i] = pi - theta[i]
  }
  sum = 0
  for (i = 0; i < nlat - 1; i++) {
    sum += wts[i]
  }
  for (i = 0; i < nlat - 1; i++) {
    wts[i] *= 2 / sum
  }
}
BEGIN {
  if (nlat == 0) nlat = 60
  rad2deg = 180 / atan2(0, -1)

  if (lgamma) {
    print gamma(1)
    print ann_gamma(1)
  }

  if (lank) {
    calc_ank(nlat, ank)
    for (i = 0; i <= nlat/2; i++) {
      printf("%5d %20.17e\n", i, ank[i])
    }
  }

  if (lcpdp) {
    cpdp(nlat, cp, dcp)
    printf("%5d %20.17e\n", 0, cp[0])
    for (i = 1; i <= (nlat+1)/2; i++) {
      printf("%5d %20.17e %20.17e\n", i, cp[i], dcp[i])
    }
  }

  if (ltpdp) {
    pi = atan2(0, -1) 
    pis2 = 0.5 * pi
    cpdp(nlat, cp, dcp)
    zero = pis2 - 0.5 * (pi / nlat)
    tpdp(nlat, zero, cp, dcp, pb, dpb)
    print pb[0], dpb[0]
  }
 
  if (lgaqd) {
    gaqd(nlat, theta, wts)
    for (i = 0; i < nlat/2; i++) {
#      printf("%5d  %20.17e %20.17e %20.17e\n", i+1, theta[i] * rad2deg, cos(theta[i]), wts[i])
      printf("%5d  %20.17e %20.17e %20.17e\n", i+1, theta[i], cos(theta[i]), wts[i])
#      printf("%5d  %20.17e  %20.17e\n", i+1, cos(theta[i]), wts[i])
    }
  }
}