MiniTensor_MatrixFunctions.h

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//                           MiniTensor Package
//
// Copyright 2016 NTESS and the MiniTensor contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER
 
#if !defined(MiniTensor_MatrixFunctions_h)
#define MiniTensor_MatrixFunctions_h
 
// Matrix functions: exp, log, sqrt families and BCH.
#include "MiniTensor_Factorizations.h"
#include "MiniTensor_Inverse.h"
#include "MiniTensor_Norms.h"
 
namespace minitensor {
 
 
template <typename T, Index N> Tensor<T, N> exp(Tensor<T, N> const &A);
 
template<typename T, Index N>
KOKKOS_INLINE_FUNCTION
Tensor<T, N>
exp_taylor(Tensor<T, N> const & A);
 
template <typename T, Index N> Tensor<T, N> exp_pade(Tensor<T, N> const &A);
 
template<typename T, Index N>
KOKKOS_INLINE_FUNCTION
Tensor<T, N>
exp_sym(Tensor<T, N> const & A);
 
template<typename T, Index N>
KOKKOS_INLINE_FUNCTION
Tensor<T, N>
exp_eig_sym(Tensor<T, N> const & A);
 
template<typename T, Index N>
KOKKOS_INLINE_FUNCTION
Tensor<T, N>
log(Tensor<T, N> const & A);
 
template<typename T, Index N>
KOKKOS_INLINE_FUNCTION
Tensor<T, N>
log_taylor(Tensor<T, N> const & A);
 
template<typename T, Index N>
KOKKOS_INLINE_FUNCTION
Tensor<T, N>
log_gregory(Tensor<T, N> const & A);
 
template<typename T, Index N>
Tensor<T, N>
log_schur(Tensor<T, N> const & A);
 
template<typename T, Index N>
KOKKOS_INLINE_FUNCTION
Tensor<T, N>
log_sym(Tensor<T, N> const & A);
 
template<typename T, Index N>
KOKKOS_INLINE_FUNCTION
Tensor<T, N>
log_eig_sym(Tensor<T, N> const & A);
 
template<typename T, Index N>
KOKKOS_INLINE_FUNCTION
Tensor<T, N>
bch(Tensor<T, N> const & v, Tensor<T, N> const & r);
 
template<typename T, Index N>
KOKKOS_INLINE_FUNCTION
Tensor<T, N>
binary_powering(Tensor<T, N> const & A, Index const exponent);
 
//
// Exponential map
//
template <typename T, Index N> Tensor<T, N> exp(Tensor<T, N> const &A) {
  return exp_pade(A);
}
 
//
// R^N exponential map by Taylor series, radius of convergence is infinity
// \param A tensor
// \return \f$ \exp A \f$
//
template<typename T, Index N>
KOKKOS_INLINE_FUNCTION
Tensor<T, N>
exp_taylor(Tensor<T, N> const & A)
{
  Index const
  max_iter = 128;
 
  T const
  tol = machine_epsilon<T>();
 
  Index const
  dimension = A.get_dimension();
 
  Tensor<T, N>
  term = identity<T, N>(dimension);
 
  // Relative error taken wrt to the first term, which is I and norm = 1
  T
  relative_error = 1.0;
 
  Tensor<T, N>
  B = term;
 
  Index
  k = 0;
 
  while (relative_error > tol && k < max_iter) {
    term = static_cast<T>(1.0 / (k + 1.0)) * term * A;
    B = B + term;
    relative_error = norm_1(term);
    ++k;
  }
 
  return B;
}
 
namespace impl {
 
//
// Scaling parameter theta for scaling and squaring exponential.
//
template<typename T>
KOKKOS_INLINE_FUNCTION
T
scaling_squaring_theta(Index const order)
{
  assert(order > 0 && order < 22);
 
  T const theta[] =
  {
      0.0e-0, 3.7e-8, 5.3e-4, 1.5e-2, 8.5e-2, 2.5e-1, 5.4e-1, 9.5e-1,
      1.5e-0, 2.1e-0, 2.8e-0, 3.6e-0, 4.5e-0, 5.4e-0, 6.3e-0, 7.3e-0,
      8.4e-0, 9.4e-0, 1.1e+1, 1.2e+1, 1.3e+1, 1.4e+1
  };
 
  return theta[order];
}
 
//
// Polynomial coefficients for Padé approximants.
//
template<typename T>
KOKKOS_INLINE_FUNCTION
T
polynomial_coefficient(Index const order, Index const index)
{
  assert(index <= order);
 
  T
  c = 0.0;
 
  switch (order) {
 
    default:
      MT_ERROR_EXIT("Wrong order in Pade' polynomial coefficient: ");
      break;
 
    case 3:
    {
      T const
      b[] = {120.0, 60.0, 12.0, 1.0};
 
      c = b[index];
    }
    break;
 
    case 5:
    {
      T const
      b[] = {30240.0, 15120.0, 3360.0, 420.0, 30.0, 1.0};
 
      c = b[index];
    }
    break;
 
    case 7:
    {
      T const
      b[] = {17297280.0, 8648640.0, 1995840.0, 277200.0, 25200.0, 1512.0,
          56.0, 1.0};
 
      c = b[index];
    }
    break;
 
    case 9:
    {
      T const
      b[] = {17643225600.0, 8821612800.0, 2075673600.0, 302702400.0,
          30270240.0, 2162160.0, 110880.0, 3960.0, 90.0, 1.0};
 
      c = b[index];
    }
    break;
 
    case 13:
    {
      T const
      b[] = {64764752532480000.0, 32382376266240000.0, 7771770303897600.0,
          1187353796428800.0, 129060195264000.0, 10559470521600.0,
          670442572800.0, 33522128640.0, 1323241920.0, 40840800.0,
          960960.0, 16380.0, 182.0, 1.0};
 
      c = b[index];
    }
    break;
 
  }
 
  return c;
}
 
//
// Padé approximant polynomial odd and even terms.
//
template <typename T, Index N>
std::pair<Tensor<T, N>, Tensor<T, N>>
pade_polynomial_terms(Tensor<T, N> const &A, Index const order) {
  Index const
  dimension = A.get_dimension();
 
  Tensor<T, N>
  B = identity<T, N>(dimension);
 
  Tensor<T, N>
  U = polynomial_coefficient<Real>(order, 1) * B;
 
  Tensor<T, N>
  V = polynomial_coefficient<Real>(order, 0) * B;
 
  Tensor<T, N> const
  A2 = A * A;
 
  for (Index i = 3; i <= order; i += 2) {
 
    B = B * A2;
 
    Tensor<T, N> const
    O = polynomial_coefficient<Real>(order, i) * B;
 
    Tensor<T, N> const
    E = polynomial_coefficient<Real>(order, i - 1) * B;
 
    U += O;
 
    V += E;
 
  }
 
  U = A * U;
 
  return std::make_pair(U, V);
}
 
//
// Coefficients for Padé approximants.
//
template<typename T>
KOKKOS_INLINE_FUNCTION
T
pade_coefficients(Index const n)
{
  T const c[] = {
    1.100343044625278e-05, 1.818617533662554e-03, 1.620628479501567e-02, 5.387353263138127e-02, 1.135280226762866e-01,
    1.866286061354130e-01, 2.642960831111435e-01, 3.402172331985299e-01, 4.108235000556820e-01, 4.745521256007768e-01,
    5.310667521178455e-01, 5.806887133441684e-01, 6.240414344012918e-01, 6.618482563071411e-01, 6.948266172489354e-01,
    7.236382701437292e-01, 7.488702930926310e-01, 7.710320825151814e-01, 7.905600074925671e-01, 8.078252198050853e-01,
    8.231422814010787e-01, 8.367774696147783e-01, 8.489562661576765e-01, 8.598698723737197e-01, 8.696807597657327e-01,
    8.785273397512191e-01, 8.865278635527148e-01, 8.937836659824918e-01, 9.003818585631236e-01, 9.063975647545747e-01,
    9.118957765024351e-01, 9.169328985287867e-01, 9.215580354375991e-01, 9.258140669835052e-01, 9.297385486977516e-01,
    9.333644683151422e-01, 9.367208829050256e-01, 9.398334570841484e-01, 9.427249190039424e-01, 9.454154478075423e-01,
    9.479230038146050e-01, 9.502636107090112e-01, 9.524515973891873e-01, 9.544998058228285e-01, 9.564197701703862e-01,
    9.582218715590143e-01, 9.599154721638511e-01, 9.615090316568806e-01, 9.630102085912245e-01, 9.644259488813590e-01,
    9.657625632018019e-01, 9.670257948457799e-01, 9.682208793510226e-01, 9.693525970039069e-01, 9.704253191689650e-01,
    9.714430492527785e-01, 9.724094589950460e-01, 9.733279206814576e-01, 9.742015357899175e-01, 9.750331605111618e-01,
    9.758254285248543e-01, 9.765807713611383e-01, 9.773014366339591e-01, 9.779895043950849e-01};
  return c[n];
}
 
template<typename T>
KOKKOS_INLINE_FUNCTION
T
gauss_legendre_abscissae(Index const m, Index const n)
{
  T x = 0.0;
  switch (m) {
    default: break;
    case 1: {
      T const c[] = {0.0};
      x           = c[n];
    } break;
    case 2: {
      T const c[] = {-0.5773502691896257645, 0.5773502691896257645};
      x           = c[n];
    } break;
 
    case 3: {
      T const c[] = {0, -0.774596669241483377, 0.774596669241483377};
      x           = c[n];
    } break;
 
    case 4: {
      T const c[] = {-0.3399810435848562648, 0.3399810435848562648, -0.8611363115940525752, 0.8611363115940525752};
      x           = c[n];
    } break;
 
    case 5: {
      T const c[] = {0, -0.5384693101056830910, 0.5384693101056830910, -0.9061798459386639927, 0.9061798459386639927};
      x           = c[n];
    } break;
 
    case 6: {
      T const c[] = {
          0.6612093864662645136613995950199,
          -0.6612093864662645136613995950199,
          -0.2386191860831969086305017216807,
          0.2386191860831969086305017216807,
          -0.9324695142031520278123015544939,
          0.9324695142031520278123015544939};
      x = c[n];
    } break;
 
    case 7: {
      T const c[] = {
          0,
          0.4058451513773971669066064120769,
          -0.4058451513773971669066064120769,
          -0.7415311855993944398638647732807,
          0.7415311855993944398638647732807,
          -0.949107912342758524526189684047851,
          0.949107912342758524526189684047851};
      x = c[n];
    } break;
 
    case 8: {
      T const c[] = {
          -0.1834346424956498049394761423601,
          0.1834346424956498049394761423601,
          -0.5255324099163289858177390491892,
          0.5255324099163289858177390491892,
          -0.7966664774136267395915539364758,
          0.7966664774136267395915539364758,
          -0.960289856497536231683560868569472,
          0.960289856497536231683560868569472};
      x = c[n];
    } break;
 
    case 9: {
      T const c[] = {
          0,
          -0.8360311073266357942994297880697,
          0.8360311073266357942994297880697,
          -0.9681602395076260898355762029036,
          0.9681602395076260898355762029036,
          -0.3242534234038089290385380146433,
          0.3242534234038089290385380146433,
          -0.6133714327005903973087020393414,
          0.6133714327005903973087020393414};
      x = c[n];
    } break;
 
    case 10: {
      T const c[] = {
          -0.1488743389816312108,
          0.1488743389816312108,
          -0.4333953941292471907,
          0.4333953941292471907,
          -0.6794095682990244062,
          0.6794095682990244062,
          -0.8650633666889845107,
          0.8650633666889845107,
          -0.9739065285171717200,
          0.9739065285171717200};
      x = c[n];
    } break;
 
    case 11: {
      T const c[] = {
          0,
          -0.2695431559523449723,
          0.2695431559523449723,
          -0.5190961292068118159,
          0.5190961292068118159,
          -0.7301520055740493240,
          0.7301520055740493240,
          -0.8870625997680952990,
          0.8870625997680952990,
          -0.978228658146056992,
          0.978228658146056992};
      x = c[n];
    } break;
 
    case 12: {
      T const c[] = {
          -0.125233408511468915,
          0.125233408511468915,
          -0.3678314989981801937,
          0.3678314989981801937,
          -0.5873179542866174472,
          0.5873179542866174472,
          -0.7699026741943046870,
          0.7699026741943046870,
          -0.9041172563704748566,
          0.9041172563704748566,
          -0.981560634246719250,
          0.981560634246719250};
      x = c[n];
    } break;
 
    case 13: {
      T const c[] = {
          0,
          -0.2304583159551347940,
          0.2304583159551347940,
          -0.4484927510364468528,
          0.4484927510364468528,
          -0.6423493394403402206,
          0.6423493394403402206,
          -0.8015780907333099127,
          0.8015780907333099127,
          -0.917598399222977965,
          0.917598399222977965,
          -0.9841830547185881494,
          0.9841830547185881494};
      x = c[n];
    } break;
 
    case 14: {
      T const c[] = {
          -0.1080549487073436620,
          0.1080549487073436620,
          -0.3191123689278897604,
          0.3191123689278897604,
          -0.5152486363581540919,
          0.5152486363581540919,
          -0.6872929048116854701,
          0.6872929048116854701,
          -0.8272013150697649931,
          0.8272013150697649931,
          -0.9284348836635735173,
          0.9284348836635735173,
          -0.9862838086968123388,
          0.9862838086968123388};
      x = c[n];
    } break;
 
    case 15: {
      T const c[] = {
          0,
          -0.2011940939974345223,
          0.2011940939974345223,
          -0.3941513470775633698,
          0.3941513470775633698,
          -0.5709721726085388475,
          0.5709721726085388475,
          -0.724417731360170047,
          0.724417731360170047,
          -0.8482065834104272162,
          0.8482065834104272162,
          -0.9372733924007059043,
          0.9372733924007059043,
          -0.9879925180204854284,
          0.9879925180204854284};
      x = c[n];
    } break;
 
    case 16: {
      T const c[] = {
          -0.0950125098376374401,
          0.0950125098376374401,
          -0.2816035507792589132,
          0.2816035507792589132,
          -0.4580167776572273863,
          0.4580167776572273863,
          -0.6178762444026437484,
          0.6178762444026437484,
          -0.7554044083550030338,
          0.7554044083550030338,
          -0.8656312023878317438,
          0.8656312023878317438,
          -0.9445750230732325760,
          0.9445750230732325760,
          -0.9894009349916499325,
          0.9894009349916499325};
      x = c[n];
    } break;
 
    case 17: {
      T const c[] = {
          0,
          -0.1784841814958478558,
          0.1784841814958478558,
          -0.3512317634538763152,
          0.3512317634538763152,
          -0.5126905370864769678,
          0.5126905370864769678,
          -0.6576711592166907658,
          0.6576711592166907658,
          -0.7815140038968014069,
          0.7815140038968014069,
          -0.8802391537269859021,
          0.8802391537269859021,
          -0.9506755217687677612,
          0.9506755217687677612,
          -0.9905754753144173356,
          0.9905754753144173356};
      x = c[n];
    } break;
 
    case 18: {
      T const c[] = {
          -0.0847750130417353012,
          0.0847750130417353012,
          -0.2518862256915055095,
          0.2518862256915055095,
          -0.4117511614628426460,
          0.4117511614628426460,
          -0.5597708310739475346,
          0.5597708310739475346,
          -0.6916870430603532078,
          0.6916870430603532078,
          -0.8037049589725231156,
          0.8037049589725231156,
          -0.8926024664975557392,
          0.8926024664975557392,
          -0.9558239495713977551,
          0.9558239495713977551,
          -0.9915651684209309467,
          0.9915651684209309467};
      x = c[n];
    } break;
 
    case 19: {
      T const c[] = {
          0,
          -0.1603586456402253758,
          0.1603586456402253758,
          -0.3165640999636298319,
          0.3165640999636298319,
          -0.4645707413759609457,
          0.4645707413759609457,
          -0.6005453046616810234,
          0.6005453046616810234,
          -0.7209661773352293786,
          0.7209661773352293786,
          -0.8227146565371428249,
          0.8227146565371428249,
          -0.9031559036148179016,
          0.9031559036148179016,
          -0.9602081521348300308,
          0.9602081521348300308,
          -0.9924068438435844031,
          0.9924068438435844031};
      x = c[n];
    } break;
 
    case 20: {
      T const c[] = {-0.0765265211334973337, 0.0765265211334973337, -0.2277858511416450780, 0.2277858511416450780,
                     -0.3737060887154195606, 0.3737060887154195606, -0.5108670019508270980, 0.5108670019508270980,
                     -0.6360536807265150254, 0.6360536807265150254, -0.7463319064601507926, 0.7463319064601507926,
                     -0.8391169718222188233, 0.8391169718222188233, -0.9122344282513259058, 0.9122344282513259058,
                     -0.9639719272779137912, 0.9639719272779137912, -0.9931285991850949247, 0.9931285991850949247};
      x           = c[n];
    } break;
 
    case 21: {
      T const c[] = {
          0,
          -0.1455618541608950909,
          0.1455618541608950909,
          -0.2880213168024010966,
          0.2880213168024010966,
          -0.4243421202074387835,
          0.4243421202074387835,
          -0.5516188358872198070,
          0.5516188358872198070,
          -0.6671388041974123193,
          0.6671388041974123193,
          -0.7684399634756779086,
          0.7684399634756779086,
          -0.8533633645833172836,
          0.8533633645833172836,
          -0.920099334150400828,
          0.920099334150400828,
          -0.9672268385663062943,
          0.9672268385663062943,
          -0.9937521706203895002,
          0.9937521706203895002};
      x = c[n];
    } break;
 
    case 22: {
      T const c[] = {-0.0697392733197222212, 0.0697392733197222212, -0.2078604266882212854, 0.2078604266882212854,
                     -0.3419358208920842251, 0.3419358208920842251, -0.4693558379867570264, 0.4693558379867570264,
                     -0.5876404035069115929, 0.5876404035069115929, -0.6944872631866827800, 0.6944872631866827800,
                     -0.7878168059792081620, 0.7878168059792081620, -0.8658125777203001365, 0.8658125777203001365,
                     -0.9269567721871740005, 0.9269567721871740005, -0.970060497835428727,  0.970060497835428727,
                     -0.9942945854823992920, 0.9942945854823992920};
      x           = c[n];
    } break;
 
    case 23: {
      T const c[] = {
          0,
          -0.133256824298466110,
          0.133256824298466110,
          -0.2641356809703449305,
          0.2641356809703449305,
          -0.3903010380302908314,
          0.3903010380302908314,
          -0.5095014778460075496,
          0.5095014778460075496,
          -0.6196098757636461563,
          0.6196098757636461563,
          -0.7186613631319501944,
          0.7186613631319501944,
          -0.8048884016188398921,
          0.8048884016188398921,
          -0.8767523582704416673,
          0.8767523582704416673,
          -0.9329710868260161023,
          0.9329710868260161023,
          -0.9725424712181152319,
          0.9725424712181152319,
          -0.994769334997552123,
          0.994769334997552123};
      x = c[n];
    } break;
 
    case 24: {
      T const c[] = {-0.0640568928626056260, 0.0640568928626056260, -0.1911188674736163091, 0.1911188674736163091,
                     -0.3150426796961633743, 0.3150426796961633743, -0.4337935076260451384, 0.4337935076260451384,
                     -0.5454214713888395356, 0.5454214713888395356, -0.6480936519369755692, 0.6480936519369755692,
                     -0.7401241915785543642, 0.7401241915785543642, -0.8200019859739029219, 0.8200019859739029219,
                     -0.8864155270044010342, 0.8864155270044010342, -0.9382745520027327585, 0.9382745520027327585,
                     -0.9747285559713094981, 0.9747285559713094981, -0.9951872199970213601, 0.9951872199970213601};
      x           = c[n];
    } break;
 
    case 25: {
      T const c[] = {
          0,
          -0.1228646926107103963,
          0.1228646926107103963,
          -0.2438668837209884320,
          0.2438668837209884320,
          -0.3611723058093878377,
          0.3611723058093878377,
          -0.4730027314457149605,
          0.4730027314457149605,
          -0.5776629302412229677,
          0.5776629302412229677,
          -0.6735663684734683644,
          0.6735663684734683644,
          -0.759259263037357630,
          0.759259263037357630,
          -0.8334426287608340014,
          0.8334426287608340014,
          -0.8949919978782753688,
          0.8949919978782753688,
          -0.9429745712289743394,
          0.9429745712289743394,
          -0.9766639214595175114,
          0.9766639214595175114,
          -0.9955569697904980979,
          0.9955569697904980979};
      x = c[n];
    } break;
 
    case 26: {
      T const c[] = {-0.0592300934293132070, 0.0592300934293132070, -0.1768588203568901839, 0.1768588203568901839,
                     -0.2920048394859568951, 0.2920048394859568951, -0.4030517551234863064, 0.4030517551234863064,
                     -0.5084407148245057176, 0.5084407148245057176, -0.6066922930176180632, 0.6066922930176180632,
                     -0.6964272604199572648, 0.6964272604199572648, -0.7763859488206788561, 0.7763859488206788561,
                     -0.8454459427884980187, 0.8454459427884980187, -0.9026378619843070742, 0.9026378619843070742,
                     -0.9471590666617142501, 0.9471590666617142501, -0.978385445956470991,  0.978385445956470991,
                     -0.9958857011456169290, 0.9958857011456169290};
      x           = c[n];
    } break;
 
    case 27: {
      T const c[] = {
          0,
          -0.1139725856095299669,
          0.1139725856095299669,
          -0.2264593654395368588,
          0.2264593654395368588,
          -0.3359939036385088997,
          0.3359939036385088997,
          -0.441148251750026880,
          0.441148251750026880,
          -0.5405515645794568949,
          0.5405515645794568949,
          -0.6329079719464951409,
          0.6329079719464951409,
          -0.7170134737394236992,
          0.7170134737394236992,
          -0.7917716390705082271,
          0.7917716390705082271,
          -0.8562079080182944903,
          0.8562079080182944903,
          -0.9094823206774911043,
          0.9094823206774911043,
          -0.9509005578147050068,
          0.9509005578147050068,
          -0.9799234759615012228,
          0.9799234759615012228,
          -0.9961792628889885669,
          0.9961792628889885669};
      x = c[n];
    } break;
 
    case 28: {
      T const c[] = {-0.0550792898840342704, 0.0550792898840342704, -0.1645692821333807712, 0.1645692821333807712,
                     -0.2720616276351780776, 0.2720616276351780776, -0.3762515160890787102, 0.3762515160890787102,
                     -0.4758742249551182610, 0.4758742249551182610, -0.5697204718114017193, 0.5697204718114017193,
                     -0.6566510940388649612, 0.6566510940388649612, -0.7356108780136317720, 0.7356108780136317720,
                     -0.8056413709171791714, 0.8056413709171791714, -0.8658925225743950489, 0.8658925225743950489,
                     -0.9156330263921320738, 0.9156330263921320738, -0.9542592806289381972, 0.9542592806289381972,
                     -0.9813031653708727536, 0.9813031653708727536, -0.996442497573954449,  0.996442497573954449};
      x           = c[n];
    } break;
 
    case 29: {
      T const c[] = {
          0,
          -0.1062782301326792301,
          0.1062782301326792301,
          -0.211352286166001074,
          0.211352286166001074,
          -0.3140316378676399349,
          0.3140316378676399349,
          -0.4131528881740086638,
          0.4131528881740086638,
          -0.5075929551242276421,
          0.5075929551242276421,
          -0.596281797138227820,
          0.596281797138227820,
          -0.6782145376026865151,
          0.6782145376026865151,
          -0.7524628517344771339,
          0.7524628517344771339,
          -0.8181854876152524449,
          0.8181854876152524449,
          -0.8746378049201027904,
          0.8746378049201027904,
          -0.9211802329530587850,
          0.9211802329530587850,
          -0.9572855957780877257,
          0.9572855957780877257,
          -0.9825455052614131748,
          0.9825455052614131748,
          -0.9966794422605965861,
          0.9966794422605965861};
      x = c[n];
    } break;
 
    case 30: {
      T const c[] = {-0.0514718425553176958, 0.0514718425553176958, -0.1538699136085835469, 0.1538699136085835469,
                     -0.2546369261678898464, 0.2546369261678898464, -0.3527047255308781134, 0.3527047255308781134,
                     -0.4470337695380891767, 0.4470337695380891767, -0.5366241481420198992, 0.5366241481420198992,
                     -0.6205261829892428611, 0.6205261829892428611, -0.697850494793315796,  0.697850494793315796,
                     -0.7677774321048261949, 0.7677774321048261949, -0.8295657623827683974, 0.8295657623827683974,
                     -0.8825605357920526815, 0.8825605357920526815, -0.926200047429274325,  0.926200047429274325,
                     -0.960021864968307512,  0.960021864968307512,  -0.9836681232797472099, 0.9836681232797472099,
                     -0.9968934840746495402, 0.9968934840746495402};
      x           = c[n];
    } break;
 
    case 31: {
      T const c[] = {
          0,
          -0.0995553121523415203,
          0.0995553121523415203,
          -0.1981211993355706287,
          0.1981211993355706287,
          -0.294718069981701616,
          0.294718069981701616,
          -0.3883859016082329430,
          0.3883859016082329430,
          -0.4781937820449024804,
          0.4781937820449024804,
          -0.563249161407149262,
          0.563249161407149262,
          -0.6427067229242603461,
          0.6427067229242603461,
          -0.7157767845868532839,
          0.7157767845868532839,
          -0.781733148416624940,
          0.781733148416624940,
          -0.8399203201462673400,
          0.8399203201462673400,
          -0.8897600299482710433,
          0.8897600299482710433,
          -0.9307569978966481649,
          0.9307569978966481649,
          -0.9625039250929496617,
          0.9625039250929496617,
          -0.9846859096651524840,
          0.9846859096651524840,
          -0.9970874818194770740,
          0.9970874818194770740};
      x = c[n];
    } break;
 
    case 32: {
      T const c[] = {-0.0483076656877383162, 0.0483076656877383162, -0.1444719615827964934, 0.1444719615827964934,
                     -0.2392873622521370745, 0.2392873622521370745, -0.3318686022821276497, 0.3318686022821276497,
                     -0.4213512761306353453, 0.4213512761306353453, -0.5068999089322293900, 0.5068999089322293900,
                     -0.5877157572407623290, 0.5877157572407623290, -0.6630442669302152009, 0.6630442669302152009,
                     -0.732182118740289680,  0.732182118740289680,  -0.7944837959679424069, 0.7944837959679424069,
                     -0.849367613732569970,  0.849367613732569970,  -0.8963211557660521239, 0.8963211557660521239,
                     -0.9349060759377396891, 0.9349060759377396891, -0.9647622555875064307, 0.9647622555875064307,
                     -0.9856115115452683354, 0.9856115115452683354, -0.9972638618494815635, 0.9972638618494815635};
      x           = c[n];
    } break;
 
    case 33: {
      T const c[] = {
          0,
          -0.0936310658547333856,
          0.0936310658547333856,
          -0.1864392988279915723,
          0.1864392988279915723,
          -0.2776090971524970294,
          0.2776090971524970294,
          -0.3663392577480733410,
          0.3663392577480733410,
          -0.4518500172724506957,
          0.4518500172724506957,
          -0.5333899047863476435,
          0.5333899047863476435,
          -0.6102423458363790273,
          0.6102423458363790273,
          -0.6817319599697427862,
          0.6817319599697427862,
          -0.7472304964495621578,
          0.7472304964495621578,
          -0.8061623562741665897,
          0.8061623562741665897,
          -0.8580096526765040646,
          0.8580096526765040646,
          -0.9023167677434335830,
          0.9023167677434335830,
          -0.9386943726111683503,
          0.9386943726111683503,
          -0.9668229096899927689,
          0.9668229096899927689,
          -0.9864557262306424881,
          0.9864557262306424881,
          -0.9974246942464552172,
          0.9974246942464552172};
      x = c[n];
    } break;
 
    case 34: {
      T const c[] = {-0.0455098219531025427, 0.0455098219531025427, -0.1361523572591829758, 0.1361523572591829758,
                     -0.2256666916164494838, 0.2256666916164494838, -0.3133110813394632474, 0.3133110813394632474,
                     -0.3983592777586459406, 0.3983592777586459406, -0.4801065451903270341, 0.4801065451903270341,
                     -0.5578755006697466427, 0.5578755006697466427, -0.6310217270805285453, 0.6310217270805285453,
                     -0.6989391132162629079, 0.6989391132162629079, -0.7610648766298730141, 0.7610648766298730141,
                     -0.8168842279009336645, 0.8168842279009336645, -0.8659346383345644692, 0.8659346383345644692,
                     -0.9078096777183244688, 0.9078096777183244688, -0.9421623974051070916, 0.9421623974051070916,
                     -0.9687082625333442817, 0.9687082625333442817, -0.9872278164063094850, 0.9872278164063094850,
                     -0.9975717537908419192, 0.9975717537908419192};
      x           = c[n];
    } break;
 
    case 35: {
      T const c[] = {
          0,
          -0.0883713432756592636,
          0.0883713432756592636,
          -0.1760510611659895699,
          0.1760510611659895699,
          -0.2623529412092960579,
          0.2623529412092960579,
          -0.3466015544308139458,
          0.3466015544308139458,
          -0.4281375415178142541,
          0.4281375415178142541,
          -0.5063227732414886150,
          0.5063227732414886150,
          -0.5805453447497645099,
          0.5805453447497645099,
          -0.650224364665890388,
          0.650224364665890388,
          -0.7148145015566287832,
          0.7148145015566287832,
          -0.7738102522869125552,
          0.7738102522869125552,
          -0.8267498990922254068,
          0.8267498990922254068,
          -0.8732191250252223315,
          0.8732191250252223315,
          -0.9128542613593176144,
          0.9128542613593176144,
          -0.9453451482078273295,
          0.9453451482078273295,
          -0.9704376160392298332,
          0.9704376160392298332,
          -0.9879357644438514980,
          0.9879357644438514980,
          -0.997706569099600297,
          0.997706569099600297};
      x = c[n];
    } break;
 
    case 36: {
      T const c[] = {-0.0430181984737086072, 0.0430181984737086072, -0.128736103809384788,  0.128736103809384788,
                     -0.2135008923168655789, 0.2135008923168655789, -0.2966849953440282705, 0.2966849953440282705,
                     -0.3776725471196892163, 0.3776725471196892163, -0.455863944433420267,  0.455863944433420267,
                     -0.5306802859262451616, 0.5306802859262451616, -0.6015676581359805350, 0.6015676581359805350,
                     -0.6680012365855210620, 0.6680012365855210620, -0.7294891715935565820, 0.7294891715935565820,
                     -0.785576230132206512,  0.785576230132206512,  -0.8358471669924753064, 0.8358471669924753064,
                     -0.8799298008903971319, 0.8799298008903971319, -0.9174977745156590660, 0.9174977745156590660,
                     -0.9482729843995075452, 0.9482729843995075452, -0.972027691049697949,  0.972027691049697949,
                     -0.9885864789022122380, 0.9885864789022122380, -0.9978304624840858361, 0.9978304624840858361};
      x           = c[n];
    } break;
 
    case 37: {
      T const c[] = {
          0,
          -0.0836704089547699019,
          0.0836704089547699019,
          -0.1667539302398519769,
          0.1667539302398519769,
          -0.248667792791365758,
          0.248667792791365758,
          -0.3288374298837069994,
          0.3288374298837069994,
          -0.4067005093183261101,
          0.4067005093183261101,
          -0.4817108778032055541,
          0.4817108778032055541,
          -0.5533423918615817812,
          0.5533423918615817812,
          -0.6210926084089244831,
          0.6210926084089244831,
          -0.6844863091309593574,
          0.6844863091309593574,
          -0.7430788339819652625,
          0.7430788339819652625,
          -0.7964592005099022933,
          0.7964592005099022933,
          -0.8442529873405559679,
          0.8442529873405559679,
          -0.8861249621554860789,
          0.8861249621554860789,
          -0.9217814374124637426,
          0.9217814374124637426,
          -0.950972343262094821,
          0.950972343262094821,
          -0.9734930300564857443,
          0.9734930300564857443,
          -0.9891859632143191866,
          0.9891859632143191866,
          -0.9979445824779136489,
          0.9979445824779136489};
      x = c[n];
    } break;
 
    case 38: {
      T const c[] = {-0.0407851479045782399, 0.0407851479045782399, -0.1220840253378674198, 0.1220840253378674198,
                     -0.2025704538921167032, 0.2025704538921167032, -0.2817088097901652613, 0.2817088097901652613,
                     -0.3589724404794350132, 0.3589724404794350132, -0.4338471694323764843, 0.4338471694323764843,
                     -0.5058347179279311032, 0.5058347179279311032, -0.5744560210478070811, 0.5744560210478070811,
                     -0.639254415829681707,  0.639254415829681707,  -0.6997986803791843559, 0.6997986803791843559,
                     -0.7556859037539706807, 0.7556859037539706807, -0.8065441676053168155, 0.8065441676053168155,
                     -0.8520350219323621888, 0.8520350219323621888, -0.8918557390046322167, 0.8918557390046322167,
                     -0.9257413320485843968, 0.9257413320485843968, -0.953466330933529595,  0.953466330933529595,
                     -0.9748463285901535076, 0.9748463285901535076, -0.9897394542663855719, 0.9897394542663855719,
                     -0.9980499305356876198, 0.9980499305356876198};
      x           = c[n];
    } break;
 
    case 39: {
      T const c[] = {
          0,
          -0.0794438046087554775,
          0.0794438046087554775,
          -0.1583853399978377999,
          0.1583853399978377999,
          -0.2363255124618357673,
          0.2363255124618357673,
          -0.3127715592481859225,
          0.3127715592481859225,
          -0.3872401639715614558,
          0.3872401639715614558,
          -0.4592605123091360486,
          0.4592605123091360486,
          -0.5283772686604374738,
          0.5283772686604374738,
          -0.5941534549572779886,
          0.5941534549572779886,
          -0.6561732134320109107,
          0.6561732134320109107,
          -0.7140444358945346791,
          0.7140444358945346791,
          -0.76740124293106349,
          0.76740124293106349,
          -0.8159062974301431043,
          0.8159062974301431043,
          -0.8592529379999061539,
          0.8592529379999061539,
          -0.8971671192929928878,
          0.8971671192929928878,
          -0.929409148486738229,
          0.929409148486738229,
          -0.9557752123246522771,
          0.9557752123246522771,
          -0.9760987093334710538,
          0.9760987093334710538,
          -0.9902515368546859836,
          0.9902515368546859836,
          -0.9981473830664329060,
          0.9981473830664329060};
      x = c[n];
    } break;
 
    case 40: {
      T const c[] = {-0.0387724175060508219, 0.0387724175060508219, -0.1160840706752552084, 0.1160840706752552084,
                     -0.1926975807013710997, 0.1926975807013710997, -0.2681521850072536811, 0.2681521850072536811,
                     -0.3419940908257584730, 0.3419940908257584730, -0.4137792043716050015, 0.4137792043716050015,
                     -0.4830758016861787129, 0.4830758016861787129, -0.5494671250951282020, 0.5494671250951282020,
                     -0.6125538896679802379, 0.6125538896679802379, -0.6719566846141795483, 0.6719566846141795483,
                     -0.727318255189927103,  0.727318255189927103,  -0.7783056514265193876, 0.7783056514265193876,
                     -0.8246122308333116631, 0.8246122308333116631, -0.8659595032122595038, 0.8659595032122595038,
                     -0.9020988069688742967, 0.9020988069688742967, -0.9328128082786765333, 0.9328128082786765333,
                     -0.9579168192137916558, 0.9579168192137916558, -0.9772599499837742626, 0.9772599499837742626,
                     -0.9907262386994570064, 0.9907262386994570064, -0.998237709710559200,  0.998237709710559200};
      x           = c[n];
    } break;
 
    case 41: {
      T const c[] = {
          0,
          -0.0756232589891629969,
          0.0756232589891629969,
          -0.1508133548639921635,
          0.1508133548639921635,
          -0.2251396056334227756,
          0.2251396056334227756,
          -0.298176277341824865,
          0.298176277341824865,
          -0.3695050226404814414,
          0.3695050226404814414,
          -0.4387172770514070885,
          0.4387172770514070885,
          -0.5054165991994060327,
          0.5054165991994060327,
          -0.569220941610215869,
          0.569220941610215869,
          -0.629764839072196320,
          0.629764839072196320,
          -0.6867015020349512895,
          0.6867015020349512895,
          -0.7397048030699261810,
          0.7397048030699261810,
          -0.7884711450474093727,
          0.7884711450474093727,
          -0.8327212004013613312,
          0.8327212004013613312,
          -0.8722015116924414088,
          0.8722015116924414088,
          -0.9066859447581011729,
          0.9066859447581011729,
          -0.935976987497853825,
          0.935976987497853825,
          -0.9599068917303462260,
          0.9599068917303462260,
          -0.9783386735610833844,
          0.9783386735610833844,
          -0.9911671096990163082,
          0.9911671096990163082,
          -0.9983215885747714415,
          0.9983215885747714415};
      x = c[n];
    } break;
 
    case 42: {
      T const c[] = {-0.036948943165351775,  0.036948943165351775,  -0.1106450272085198683, 0.1106450272085198683,
                     -0.1837368065648545508, 0.1837368065648545508, -0.2558250793428790839, 0.2558250793428790839,
                     -0.3265161244654115121, 0.3265161244654115121, -0.3954238520429750576, 0.3954238520429750576,
                     -0.4621719120704219297, 0.4621719120704219297, -0.5263957499311922875, 0.5263957499311922875,
                     -0.5877445974851093228, 0.5877445974851093228, -0.6458833888692478339, 0.6458833888692478339,
                     -0.700494590556171213,  0.700494590556171213,  -0.7512799356894804895, 0.7512799356894804895,
                     -0.7979620532554874132, 0.7979620532554874132, -0.8402859832618169009, 0.8402859832618169009,
                     -0.8780205698121727427, 0.8780205698121727427, -0.9109597249041274525, 0.9109597249041274525,
                     -0.9389235573549881785, 0.9389235573549881785, -0.9617593653382044887, 0.9617593653382044887,
                     -0.979342508063748193,  0.979342508063748193,  -0.9915772883408609197, 0.9915772883408609197,
                     -0.9983996189900624150, 0.9983996189900624150};
      x           = c[n];
    } break;
 
    case 43: {
      T const c[] = {
          0,
          -0.0721529908745862354,
          0.0721529908745862354,
          -0.1439298095107133107,
          0.1439298095107133107,
          -0.2149562448605182090,
          0.2149562448605182090,
          -0.2848619980329136271,
          0.2848619980329136271,
          -0.3532826128643038066,
          0.3532826128643038066,
          -0.4198613760292692524,
          0.4198613760292692524,
          -0.4842511767857347240,
          0.4842511767857347240,
          -0.5461163166600847191,
          0.5461163166600847191,
          -0.605134259639600935,
          0.605134259639600935,
          -0.6609973137514981331,
          0.6609973137514981331,
          -0.713414235268957054,
          0.713414235268957054,
          -0.7621117471949551214,
          0.7621117471949551214,
          -0.8068359641369386352,
          0.8068359641369386352,
          -0.8473537162093150489,
          0.8473537162093150489,
          -0.8834537652186168633,
          0.8834537652186168633,
          -0.914947907206138729,
          0.914947907206138729,
          -0.9416719568476378618,
          0.9416719568476378618,
          -0.9634866130140799934,
          0.9634866130140799934,
          -0.9802782209802553315,
          0.9802782209802553315,
          -0.9919595575932441464,
          0.9919595575932441464,
          -0.9984723322425077135,
          0.9984723322425077135};
      x = c[n];
    } break;
 
    case 44: {
      T const c[] = {-0.0352892369641353590, 0.0352892369641353590, -0.105691901708653247,  0.105691901708653247,
                     -0.1755680147755167857, 0.1755680147755167857, -0.2445694569282012515, 0.2445694569282012515,
                     -0.3123524665027858122, 0.3123524665027858122, -0.3785793520147071325, 0.3785793520147071325,
                     -0.4429201745254114838, 0.4429201745254114838, -0.5050543913882023179, 0.5050543913882023179,
                     -0.5646724531854707684, 0.5646724531854707684, -0.6214773459035758478, 0.6214773459035758478,
                     -0.6751860706661223653, 0.6751860706661223653, -0.7255310536607170026, 0.7255310536607170026,
                     -0.7722614792487558990, 0.7722614792487558990, -0.8151445396451350104, 0.8151445396451350104,
                     -0.8539665950047103787, 0.8539665950047103787, -0.8885342382860432023, 0.8885342382860432023,
                     -0.91867525998417577,   0.91867525998417577,   -0.9442395091181940992, 0.9442395091181940992,
                     -0.9650996504224931393, 0.9650996504224931393, -0.9811518330779139666, 0.9811518330779139666,
                     -0.9923163921385158084, 0.9923163921385158084, -0.9985402006367742249, 0.9985402006367742249};
      x           = c[n];
    } break;
 
    case 45: {
      T const c[] = {
          0,
          -0.0689869801631441724,
          0.0689869801631441724,
          -0.1376452059832530287,
          0.1376452059832530287,
          -0.2056474897832637457,
          0.2056474897832637457,
          -0.2726697697523775606,
          0.2726697697523775606,
          -0.338392654250602161,
          0.338392654250602161,
          -0.4025029438585419140,
          0.4025029438585419140,
          -0.4646951239196350985,
          0.4646951239196350985,
          -0.5246728204629160670,
          0.5246728204629160670,
          -0.5821502125693531866,
          0.5821502125693531866,
          -0.6368533944532233592,
          0.6368533944532233592,
          -0.6885216807712005252,
          0.6885216807712005252,
          -0.7369088489454903526,
          0.7369088489454903526,
          -0.7817843125939062913,
          0.7817843125939062913,
          -0.8229342205020863370,
          0.8229342205020863370,
          -0.8601624759606642253,
          0.8601624759606642253,
          -0.8932916717532417384,
          0.8932916717532417384,
          -0.9221639367190003880,
          0.9221639367190003880,
          -0.9466416909956290617,
          0.9466416909956290617,
          -0.9666083103968946047,
          0.9666083103968946047,
          -0.9819687150345405682,
          0.9819687150345405682,
          -0.992649998447203741,
          0.992649998447203741,
          -0.9986036451819366381,
          0.9986036451819366381};
      x = c[n];
    } break;
 
    case 46: {
      T const c[] = {-0.0337721900160520415, 0.0337721900160520415, -0.101162475305584239,  0.101162475305584239,
                     -0.1680911794671035286, 0.1680911794671035286, -0.2342529222062697686, 0.2342529222062697686,
                     -0.29934582270187001,   0.29934582270187001,   -0.3630728770209957101, 0.3630728770209957101,
                     -0.4251433132828283973, 0.4251433132828283973, -0.4852739183881646627, 0.4852739183881646627,
                     -0.5431903302618026352, 0.5431903302618026352, -0.598628289712715153,  0.598628289712715153,
                     -0.6513348462019977151, 0.6513348462019977151, -0.7010695120204056975, 0.7010695120204056975,
                     -0.747605359615666054,  0.747605359615666054,  -0.7907300570752742551, 0.7907300570752742551,
                     -0.830246837066066053,  0.830246837066066053,  -0.8659753948668580629, 0.8659753948668580629,
                     -0.8977527115339419657, 0.8977527115339419657, -0.9254337988067539509, 0.9254337988067539509,
                     -0.9488923634460897956, 0.9488923634460897956, -0.9680213918539919427, 0.9680213918539919427,
                     -0.9827336698041668634, 0.9827336698041668634, -0.9929623489061743640, 0.9929623489061743640,
                     -0.9986630421338179811, 0.9986630421338179811};
      x           = c[n];
    } break;
 
    case 47: {
      T const c[] = {
          0,
          -0.0660869239163556751,
          0.0660869239163556751,
          -0.1318848665545148970,
          0.1318848665545148970,
          -0.1971061102791118079,
          0.1971061102791118079,
          -0.2614654592149745703,
          0.2614654592149745703,
          -0.3246814863377359022,
          0.3246814863377359022,
          -0.3864777640846671395,
          0.3864777640846671395,
          -0.4465840731048557027,
          0.4465840731048557027,
          -0.5047375838635779197,
          0.5047375838635779197,
          -0.5606840059346641944,
          0.5606840059346641944,
          -0.6141786999563736085,
          0.6141786999563736085,
          -0.6649877473903327291,
          0.6649877473903327291,
          -0.7128889734090643016,
          0.7128889734090643016,
          -0.7576729184454386335,
          0.7576729184454386335,
          -0.7991437541677419429,
          0.7991437541677419429,
          -0.8371201398999021212,
          0.8371201398999021212,
          -0.8714360157968963169,
          0.8714360157968963169,
          -0.9019413294385253568,
          0.9019413294385253568,
          -0.9285026930123606481,
          0.9285026930123606481,
          -0.9510039692577084425,
          0.9510039692577084425,
          -0.9693467873265644971,
          0.9693467873265644971,
          -0.9834510030716237087,
          0.9834510030716237087,
          -0.9932552109877686346,
          0.9932552109877686346,
          -0.998718728584212109,
          0.998718728584212109};
      x = c[n];
    } break;
 
    case 48: {
      T const c[] = {-0.0323801709628693620, 0.0323801709628693620, -0.097004699209462698,  0.097004699209462698,
                     -0.1612223560688917180, 0.1612223560688917180, -0.2247637903946890612, 0.2247637903946890612,
                     -0.287362487355455576,  0.287362487355455576,  -0.3487558862921607381, 0.3487558862921607381,
                     -0.4086864819907167299, 0.4086864819907167299, -0.4669029047509584045, 0.4669029047509584045,
                     -0.5231609747222330336, 0.5231609747222330336, -0.5772247260839727038, 0.5772247260839727038,
                     -0.6288673967765136239, 0.6288673967765136239, -0.677872379632663905,  0.677872379632663905,
                     -0.7240341309238146546, 0.7240341309238146546, -0.7671590325157403392, 0.7671590325157403392,
                     -0.8070662040294426270, 0.8070662040294426270, -0.8435882616243935307, 0.8435882616243935307,
                     -0.8765720202742478859, 0.8765720202742478859, -0.9058791367155696728, 0.9058791367155696728,
                     -0.9313866907065543331, 0.9313866907065543331, -0.9529877031604308607, 0.9529877031604308607,
                     -0.9705915925462472504, 0.9705915925462472504, -0.9841245837228268577, 0.9841245837228268577,
                     -0.9935301722663507575, 0.9935301722663507575, -0.9987710072524261186, 0.9987710072524261186};
      x           = c[n];
    } break;
 
    case 49: {
      T const c[] = {
          0,
          -0.0634206849826867860,
          0.0634206849826867860,
          -0.1265859972696720510,
          0.1265859972696720510,
          -0.189241592461813586,
          0.189241592461813586,
          -0.2511351786125772735,
          0.2511351786125772735,
          -0.3120175321197487622,
          0.3120175321197487622,
          -0.3716435012622848888,
          0.3716435012622848888,
          -0.4297729933415765246,
          0.4297729933415765246,
          -0.4861719414524920421,
          0.4861719414524920421,
          -0.5406132469917260665,
          0.5406132469917260665,
          -0.5928776941089007124,
          0.5928776941089007124,
          -0.642754832419237664,
          0.642754832419237664,
          -0.6900438244251321135,
          0.6900438244251321135,
          -0.7345542542374026962,
          0.7345542542374026962,
          -0.7761068943454466350,
          0.7761068943454466350,
          -0.8145344273598554315,
          0.8145344273598554315,
          -0.8496821198441657010,
          0.8496821198441657010,
          -0.8814084455730089100,
          0.8814084455730089100,
          -0.9095856558280732852,
          0.9095856558280732852,
          -0.934100294755810149,
          0.934100294755810149,
          -0.9548536586741372335,
          0.9548536586741372335,
          -0.9717622009015553801,
          0.9717622009015553801,
          -0.9847578959142130043,
          0.9847578959142130043,
          -0.9937886619441677907,
          0.9937886619441677907,
          -0.9988201506066353793,
          0.9988201506066353793};
      x = c[n];
    } break;
 
    case 50: {
      T const c[] = {-0.0310983383271888761, 0.0310983383271888761, -0.0931747015600861408, 0.0931747015600861408,
                     -0.1548905899981459020, 0.1548905899981459020, -0.2160072368760417568, 0.2160072368760417568,
                     -0.2762881937795319903, 0.2762881937795319903, -0.3355002454194373568, 0.3355002454194373568,
                     -0.3934143118975651273, 0.3934143118975651273, -0.4498063349740387891, 0.4498063349740387891,
                     -0.504458144907464201,  0.504458144907464201,  -0.5571583045146500543, 0.5571583045146500543,
                     -0.6077029271849502391, 0.6077029271849502391, -0.6558964656854393607, 0.6558964656854393607,
                     -0.701552468706822251,  0.701552468706822251,  -0.7444943022260685382, 0.7444943022260685382,
                     -0.7845558329003992639, 0.7845558329003992639, -0.8215820708593359483, 0.8215820708593359483,
                     -0.8554297694299460846, 0.8554297694299460846, -0.8859679795236130486, 0.8859679795236130486,
                     -0.9130785566557918930, 0.9130785566557918930, -0.9366566189448779337, 0.9366566189448779337,
                     -0.9566109552428079429, 0.9566109552428079429, -0.9728643851066920737, 0.9728643851066920737,
                     -0.9853540840480058823, 0.9853540840480058823, -0.9940319694320907125, 0.9940319694320907125,
                     -0.9988664044200710501, 0.9988664044200710501};
      x           = c[n];
    } break;
 
    case 51: {
      T const c[] = {
          0,
          -0.0609611001505787247,
          0.0609611001505787247,
          -0.121695421018888766,
          0.121695421018888766,
          -0.1819770269570775453,
          0.1819770269570775453,
          -0.2415816664477987038,
          0.2415816664477987038,
          -0.3002876063353319395,
          0.3002876063353319395,
          -0.3578764566884095097,
          0.3578764566884095097,
          -0.4141339832263038779,
          0.4141339832263038779,
          -0.4688509042860410636,
          0.4688509042860410636,
          -0.521823669366185842,
          0.521823669366185842,
          -0.5728552163513038365,
          0.5728552163513038365,
          -0.6217557046007232737,
          0.6217557046007232737,
          -0.6683432211753700868,
          0.6683432211753700868,
          -0.7124444575770366445,
          0.7124444575770366445,
          -0.7538953544853755257,
          0.7538953544853755257,
          -0.79254171209938120,
          0.79254171209938120,
          -0.8282397638230648328,
          0.8282397638230648328,
          -0.8608567111822923714,
          0.8608567111822923714,
          -0.8902712180295273032,
          0.8902712180295273032,
          -0.9163738623097802308,
          0.9163738623097802308,
          -0.9390675440029623834,
          0.9390675440029623834,
          -0.9582678486139081945,
          0.9582678486139081945,
          -0.9739033680193238672,
          0.9739033680193238672,
          -0.9859159917359029965,
          0.9859159917359029965,
          -0.9942612604367525746,
          0.9942612604367525746,
          -0.9989099908489034951,
          0.9989099908489034951};
      x = c[n];
    } break;
 
    case 52: {
      T const c[] = {-0.0299141097973387660, 0.0299141097973387660, -0.089635244648900565,  0.089635244648900565,
                     -0.1490355086069491804, 0.1490355086069491804, -0.2079022641563660596, 0.2079022641563660596,
                     -0.2660247836050018274, 0.2660247836050018274, -0.3231950034348078255, 0.3231950034348078255,
                     -0.3792082691160936692, 0.3792082691160936692, -0.4338640677187616703, 0.4338640677187616703,
                     -0.4869667456980960777, 0.4869667456980960777, -0.5383262092858274383, 0.5383262092858274383,
                     -0.5877586049795790699, 0.5877586049795790699, -0.6350869776952459242, 0.6350869776952459242,
                     -0.6801419042271677020, 0.6801419042271677020, -0.7227620997499831936, 0.7227620997499831936,
                     -0.7627949951937449602, 0.7627949951937449602, -0.800097283430468324,  0.800097283430468324,
                     -0.8345354323267345349, 0.8345354323267345349, -0.8659861628460675852, 0.8659861628460675852,
                     -0.8943368905344953225, 0.8943368905344953225, -0.9194861289164245398, 0.9194861289164245398,
                     -0.9413438536413590568, 0.9413438536413590568, -0.9598318269330865525, 0.9598318269330865525,
                     -0.9748838842217445031, 0.9748838842217445031, -0.9864461956515498406, 0.9864461956515498406,
                     -0.9944775909292160292, 0.9944775909292160292, -0.9989511111039502780, 0.9989511111039502780};
      x           = c[n];
    } break;
 
    case 53: {
      T const c[] = {
          0,
          -0.058685054300259465,
          0.058685054300259465,
          -0.1171678090719551501,
          0.1171678090719551501,
          -0.1752466621553257507,
          0.1752466621553257507,
          -0.232721403724272593,
          0.232721403724272593,
          -0.289393906451626206,
          0.289393906451626206,
          -0.3450688084957223566,
          0.3450688084957223566,
          -0.3995541869539529773,
          0.3995541869539529773,
          -0.4526622194618457913,
          0.4526622194618457913,
          -0.5042098316571334370,
          0.5042098316571334370,
          -0.5540193282770678810,
          0.5540193282770678810,
          -0.6019190057137693274,
          0.6019190057137693274,
          -0.6477437439165100687,
          0.6477437439165100687,
          -0.6913355756013667235,
          0.6913355756013667235,
          -0.7325442308075102537,
          0.7325442308075102537,
          -0.7712276549255323078,
          0.7712276549255323078,
          -0.807252498416895478,
          0.807252498416895478,
          -0.8404945765458013754,
          0.8404945765458013754,
          -0.8708392975582413516,
          0.8708392975582413516,
          -0.8981820578754266259,
          0.8981820578754266259,
          -0.9224286030428121282,
          0.9224286030428121282,
          -0.9434953534644418790,
          0.9434953534644418790,
          -0.9613096946231363323,
          0.9613096946231363323,
          -0.9758102337149845816,
          0.9758102337149845816,
          -0.9869470350233715217,
          0.9869470350233715217,
          -0.9946819193080070786,
          0.9946819193080070786,
          -0.9989899477763282271,
          0.9989899477763282271};
      x = c[n];
    } break;
 
    case 54: {
      T const c[] = {-0.0288167481993417776, 0.0288167481993417776, -0.086354518263248215,  0.086354518263248215,
                     -0.1436054273162561539, 0.1436054273162561539, -0.2003792936062135697, 0.2003792936062135697,
                     -0.2564875200699973000, 0.2564875200699973000, -0.3117437208344682288, 0.3117437208344682288,
                     -0.3659643403721911819, 0.3659643403721911819, -0.418969263255204528,  0.418969263255204528,
                     -0.4705824124813822836, 0.4705824124813822836, -0.5206323343859330733, 0.5206323343859330733,
                     -0.5689527681952094297, 0.5689527681952094297, -0.6153831983311273707, 0.6153831983311273707,
                     -0.6597693876319831246, 0.6597693876319831246, -0.7019638897191729193, 0.7019638897191729193,
                     -0.7418265388091843162, 0.7418265388091843162, -0.7792249153462540215, 0.7792249153462540215,
                     -0.8140347859135678354, 0.8140347859135678354, -0.8461405159707729494, 0.8461405159707729494,
                     -0.8754354540655689394, 0.8754354540655689394, -0.901822286284701580,  0.901822286284701580,
                     -0.9252133598666514862, 0.9252133598666514862, -0.9455309751649958537, 0.9455309751649958537,
                     -0.9627076457859235832, 0.9627076457859235832, -0.9766863288579032372, 0.9766863288579032372,
                     -0.9874206373973435585, 0.9874206373973435585, -0.9948751170183388849, 0.9948751170183388849,
                     -0.9990266668673409838, 0.9990266668673409838};
      x           = c[n];
    } break;
 
    case 55: {
      T const c[] = {
          0,
          -0.0565727538183367763,
          0.0565727538183367763,
          -0.1129642880593292665,
          0.1129642880593292665,
          -0.1689939636468732082,
          0.1689939636468732082,
          -0.2244823006478454834,
          0.2244823006478454834,
          -0.2792515532008065385,
          0.2792515532008065385,
          -0.3331262788900238851,
          0.3331262788900238851,
          -0.3859339007409794297,
          0.3859339007409794297,
          -0.4375052600371745918,
          0.4375052600371745918,
          -0.4876751581874740972,
          0.4876751581874740972,
          -0.5362828859083432967,
          0.5362828859083432967,
          -0.5831727380260321029,
          0.5831727380260321029,
          -0.6281945122499281400,
          0.6281945122499281400,
          -0.6712039903198263957,
          0.6712039903198263957,
          -0.7120633999866378389,
          0.7120633999866378389,
          -0.7506418563480219086,
          0.7506418563480219086,
          -0.7868157811276223658,
          0.7868157811276223658,
          -0.8204692985593209124,
          0.8204692985593209124,
          -0.8514946066171544714,
          0.8514946066171544714,
          -0.8797923224198955060,
          0.8797923224198955060,
          -0.905271800744000025,
          0.905271800744000025,
          -0.9278514247207916968,
          0.9278514247207916968,
          -0.9474588680412107418,
          0.9474588680412107418,
          -0.9640313285931351987,
          0.9640313285931351987,
          -0.9775157355039892088,
          0.9775157355039892088,
          -0.9878689411988891985,
          0.9878689411988891985,
          -0.9950579778474118750,
          0.9950579778474118750,
          -0.9990614195648185414,
          0.9990614195648185414};
      x = c[n];
    } break;
 
    case 56: {
      T const c[] = {-0.0277970352872754370, 0.0277970352872754370, -0.0833051868224353744, 0.0833051868224353744,
                     -0.1385558468103762420, 0.1385558468103762420, -0.1933782386352752582, 0.1933782386352752582,
                     -0.2476029094343372039, 0.2476029094343372039, -0.3010622538672206690, 0.3010622538672206690,
                     -0.3535910321749545209, 0.3535910321749545209, -0.405026880927091278,  0.405026880927091278,
                     -0.455210814878459578,  0.455210814878459578,  -0.5039877183843817141, 0.5039877183843817141,
                     -0.5512068248555346187, 0.5512068248555346187, -0.5967221827706633201, 0.5967221827706633201,
                     -0.6403931068070068942, 0.6403931068070068942, -0.6820846126944704555, 0.6820846126944704555,
                     -0.7216678344501880835, 0.7216678344501880835, -0.7590204227051289022, 0.7590204227051289022,
                     -0.7940269228938664980, 0.7940269228938664980, -0.8265791321428816516, 0.8265791321428816516,
                     -0.8565764337627486354, 0.8565764337627486354, -0.8839261083278275407, 0.8839261083278275407,
                     -0.9085436204206554908, 0.9085436204206554908, -0.9303528802474963005, 0.9303528802474963005,
                     -0.9492864795619626356, 0.9492864795619626356, -0.9652859019054901836, 0.9652859019054901836,
                     -0.9783017091402563833, 0.9783017091402563833, -0.9882937155401615110, 0.9882937155401615110,
                     -0.9952312260810697472, 0.9952312260810697472, -0.999094343801465584,  0.999094343801465584};
      x           = c[n];
    } break;
 
    case 57: {
      T const c[] = {
          0,
          -0.0546071510016468242,
          0.0546071510016468242,
          -0.1090513328087878009,
          0.1090513328087878009,
          -0.1631700625912642510,
          0.1631700625912642510,
          -0.2168018287961240364,
          0.2168018287961240364,
          -0.2697865731618387657,
          0.2697865731618387657,
          -0.3219661683953786405,
          0.3219661683953786405,
          -0.3731848900865944585,
          0.3731848900865944585,
          -0.4232898814515639509,
          0.4232898814515639509,
          -0.4721316095179757095,
          0.4721316095179757095,
          -0.5195643113911876063,
          0.5195643113911876063,
          -0.5654464292692367590,
          0.5654464292692367590,
          -0.6096410329087153654,
          0.6096410329087153654,
          -0.6520162282809768912,
          0.6520162282809768912,
          -0.692445551199517739,
          0.692445551199517739,
          -0.7308083447445233228,
          0.7308083447445233228,
          -0.7669901193594501954,
          0.7669901193594501954,
          -0.8008828945472182420,
          0.8008828945472182420,
          -0.8323855211504391208,
          0.8323855211504391208,
          -0.8614039832620469447,
          0.8614039832620469447,
          -0.8878516788822213295,
          0.8878516788822213295,
          -0.9116496785213912127,
          0.9116496785213912127,
          -0.9327269610671016961,
          0.9327269610671016961,
          -0.9510206264478767419,
          0.9510206264478767419,
          -0.9664760851718866791,
          0.9664760851718866791,
          -0.9790472267094687137,
          0.9790472267094687137,
          -0.9886965776502220488,
          0.9886965776502220488,
          -0.9953955236784303111,
          0.9953955236784303111,
          -0.9991255656252628505,
          0.9991255656252628505};
      x = c[n];
    } break;
 
    case 58: {
      T const c[] = {-0.0268470123659423558, 0.0268470123659423558, -0.0804636302141427293, 0.0804636302141427293,
                     -0.1338482505954668570, 0.1338482505954668570, -0.1868469518357613213, 0.1868469518357613213,
                     -0.2393069249661534544, 0.2393069249661534544, -0.2910769143111091895, 0.2910769143111091895,
                     -0.3420076535979952612, 0.3420076535979952612, -0.3919522963307531503, 0.3919522963307531503,
                     -0.4407668391868395651, 0.4407668391868395651, -0.4883105372167184636, 0.4883105372167184636,
                     -0.5344463096488475863, 0.5344463096488475863, -0.5790411351302250304, 0.5790411351302250304,
                     -0.6219664352630791110, 0.6219664352630791110, -0.6630984453321252664, 0.6630984453321252664,
                     -0.7023185711539081134, 0.7023185711539081134, -0.7395137310200422678, 0.7395137310200422678,
                     -0.7745766817496527452, 0.7745766817496527452, -0.8074063279130881410, 0.8074063279130881410,
                     -0.8379080133393733163, 0.8379080133393733163, -0.8659937940748074792, 0.8659937940748074792,
                     -0.8915826920220301763, 0.8915826920220301763, -0.9146009285643525406, 0.9146009285643525406,
                     -0.9349821375882593484, 0.9349821375882593484, -0.9526675575188690914, 0.9526675575188690914,
                     -0.9676062025029240901, 0.9676062025029240901, -0.9797550146943503091, 0.9797550146943503091,
                     -0.9890790082484426364, 0.9890790082484426364, -0.9955514765972909026, 0.9955514765972909026,
                     -0.999155200407386606,  0.999155200407386606};
      x           = c[n];
    } break;
 
    case 59: {
      T const c[] = {
          0,
          -0.0527734840883100039,
          0.0527734840883100039,
          -0.1053998790163441438,
          0.1053998790163441438,
          -0.1577325055878579681,
          0.1577325055878579681,
          -0.2096255033920365449,
          0.2096255033920365449,
          -0.260934237342811711,
          0.260934237342811711,
          -0.3115157008030137003,
          0.3115157008030137003,
          -0.3612289141697948099,
          0.3612289141697948099,
          -0.4099353178104189667,
          0.4099353178104189667,
          -0.4574991582532666902,
          0.4574991582532666902,
          -0.5037878665577179787,
          0.5037878665577179787,
          -0.5486724278083963843,
          0.5486724278083963843,
          -0.5920277407040301444,
          0.5920277407040301444,
          -0.6337329662388500975,
          0.6337329662388500975,
          -0.6736718645049372270,
          0.6736718645049372270,
          -0.7117331186771977315,
          0.7117331186771977315,
          -0.747810645278640231,
          0.747810645278640231,
          -0.7818038898623609056,
          0.7818038898623609056,
          -0.8136181072882115714,
          0.8136181072882115714,
          -0.8431646258168722014,
          0.8431646258168722014,
          -0.8703610942928822609,
          0.8703610942928822609,
          -0.895131711743472085,
          0.895131711743472085,
          -0.9174074387881552813,
          0.9174074387881552813,
          -0.9371261903534538594,
          0.9371261903534538594,
          -0.9542330093769510558,
          0.9542330093769510558,
          -0.9686802216817815313,
          0.9686802216817815313,
          -0.9804275739567156884,
          0.9804275739567156884,
          -0.9894423651337309317,
          0.9894423651337309317,
          -0.995699640383245964,
          0.995699640383245964,
          -0.999183353909294683,
          0.999183353909294683};
      x = c[n];
    } break;
 
    case 60: {
      T const c[] = {-0.0259597723012477985, 0.0259597723012477985, -0.0778093339495365694, 0.0778093339495365694,
                     -0.1294491353969450031, 0.1294491353969450031, -0.1807399648734254172, 0.1807399648734254172,
                     -0.2315435513760293380, 0.2315435513760293380, -0.2817229374232616916, 0.2817229374232616916,
                     -0.3311428482684481942, 0.3311428482684481942, -0.3796700565767979771, 0.3796700565767979771,
                     -0.4271737415830783893, 0.4271737415830783893, -0.4735258417617071111, 0.4735258417617071111,
                     -0.5186014000585697474, 0.5186014000585697474, -0.562278900753944539,  0.562278900753944539,
                     -0.6044405970485103634, 0.6044405970485103634, -0.6449728284894770678, 0.6449728284894770678,
                     -0.6837663273813554372, 0.6837663273813554372, -0.7207165133557303994, 0.7207165133557303994,
                     -0.7557237753065856868, 0.7557237753065856868, -0.7886937399322640545, 0.7886937399322640545,
                     -0.8195375261621457593, 0.8195375261621457593, -0.8481719847859296324, 0.8481719847859296324,
                     -0.8745199226468983151, 0.8745199226468983151, -0.8985103108100459419, 0.8985103108100459419,
                     -0.9200784761776275528, 0.9200784761776275528, -0.9391662761164232494, 0.9391662761164232494,
                     -0.9557222558399961073, 0.9557222558399961073, -0.9697017887650527337, 0.9697017887650527337,
                     -0.9810672017525981856, 0.9810672017525981856, -0.9897878952222217173, 0.9897878952222217173,
                     -0.9958405251188381738, 0.9958405251188381738, -0.9992101232274360220, 0.9992101232274360220};
      x           = c[n];
    } break;
 
    case 61: {
      T const c[] = {
          0,
          -0.0510589067079743493,
          0.0510589067079743493,
          -0.1019846065622740689,
          0.1019846065622740689,
          -0.1526442402308153005,
          0.1526442402308153005,
          -0.2029056425180584992,
          0.2029056425180584992,
          -0.2526376871690534958,
          0.2526376871690534958,
          -0.3017106289630307126,
          0.3017106289630307126,
          -0.3499964422040668345,
          0.3499964422040668345,
          -0.3973691547257566091,
          0.3973691547257566091,
          -0.4437051765385316019,
          0.4437051765385316019,
          -0.4888836222622521188,
          0.4888836222622521188,
          -0.5327866265029252656,
          0.5327866265029252656,
          -0.5752996513508306186,
          0.5752996513508306186,
          -0.6163117851979217247,
          0.6163117851979217247,
          -0.6557160320950708716,
          0.6557160320950708716,
          -0.6934095908944911554,
          0.6934095908944911554,
          -0.7292941234494651096,
          0.7292941234494651096,
          -0.7632760111723121971,
          0.7632760111723121971,
          -0.7952665992823596491,
          0.7952665992823596491,
          -0.8251824281086599506,
          0.8251824281086599506,
          -0.8529454508476634455,
          0.8529454508476634455,
          -0.8784832372148810324,
          0.8784832372148810324,
          -0.901729162474001170,
          0.901729162474001170,
          -0.9226225813829552612,
          0.9226225813829552612,
          -0.9411089866813611474,
          0.9411089866813611474,
          -0.9571401519129840913,
          0.9571401519129840913,
          -0.9706742588331829082,
          0.9706742588331829082,
          -0.9816760112840370796,
          0.9816760112840370796,
          -0.9901167452325170509,
          0.9901167452325170509,
          -0.9959745998151202342,
          0.9959745998151202342,
          -0.9992355976313634717,
          0.9992355976313634717};
      x = c[n];
    } break;
 
    case 62: {
      T const c[] = {-0.0251292914218206147, 0.0251292914218206147, -0.0753243954962343327, 0.0753243954962343327,
                     -0.1253292236158968086, 0.1253292236158968086, -0.175017459249015628,  0.175017459249015628,
                     -0.2242635856041655316, 0.2242635856041655316, -0.2729432026967263431, 0.2729432026967263431,
                     -0.3209333415941940040, 0.3209333415941940040, -0.3681127750465645296, 0.3681127750465645296,
                     -0.4143623237171260481, 0.4143623237171260481, -0.4595651572401133952, 0.4595651572401133952,
                     -0.5036070893447559559, 0.5036070893447559559, -0.5463768663002510958, 0.5463768663002510958,
                     -0.5877664479530873380, 0.5877664479530873380, -0.6276712806468851807, 0.6276712806468851807,
                     -0.665990561335479446,  0.665990561335479446,  -0.7026274922222970551, 0.7026274922222970551,
                     -0.7374895252831567498, 0.7374895252831567498, -0.7704885960554193189, 0.7704885960554193189,
                     -0.8015413461039763715, 0.8015413461039763715, -0.8305693336040048513, 0.8305693336040048513,
                     -0.8574992315120709228, 0.8574992315120709228, -0.882263012831897363,  0.882263012831897363,
                     -0.9047981225210934657, 0.9047981225210934657, -0.9250476356362037552, 0.9250476356362037552,
                     -0.9429604013923285038, 0.9429604013923285038, -0.95849117297392709,   0.95849117297392709,
                     -0.9716007233716518064, 0.9716007233716518064, -0.9822559490972366494, 0.9822559490972366494,
                     -0.9904299711892903524, 0.9904299711892903524, -0.9961022963162671328, 0.9961022963162671328,
                     -0.9992598593087770296, 0.9992598593087770296};
      x           = c[n];
    } break;
 
    case 63: {
      T const c[] = {
          0,
          -0.0494521871161596272,
          0.0494521871161596272,
          -0.0987833564469452795,
          0.0987833564469452795,
          -0.1478727863578719685,
          0.1478727863578719685,
          -0.1966003467915066845,
          0.1966003467915066845,
          -0.2448467932459533627,
          0.2448467932459533627,
          -0.2924940585862514400,
          0.2924940585862514400,
          -0.3394255419745844024,
          0.3394255419745844024,
          -0.385526394212247892,
          0.385526394212247892,
          -0.4306837987951116006,
          0.4306837987951116006,
          -0.4747872479948043999,
          0.4747872479948043999,
          -0.5177288132900332481,
          0.5177288132900332481,
          -0.5594034094862850132,
          0.5594034094862850132,
          -0.5997090518776252357,
          0.5997090518776252357,
          -0.6385471058213653850,
          0.6385471058213653850,
          -0.6758225281149860901,
          0.6758225281149860901,
          -0.7114440995848458078,
          0.7114440995848458078,
          -0.7453246483178474178,
          0.7453246483178474178,
          -0.7773812629903723355,
          0.7773812629903723355,
          -0.8075354957734567600,
          0.8075354957734567600,
          -0.8357135543195028434,
          0.8357135543195028434,
          -0.8618464823641237195,
          0.8618464823641237195,
          -0.8858703285078534262,
          0.8858703285078534262,
          -0.907726302778531558,
          0.907726302778531558,
          -0.9273609206218432054,
          0.9273609206218432054,
          -0.9447261340410098029,
          0.9447261340410098029,
          -0.959779449758941927,
          0.959779449758941927,
          -0.9724840346975700228,
          0.9724840346975700228,
          -0.9828088105937272348,
          0.9828088105937272348,
          -0.9907285468921894668,
          0.9907285468921894668,
          -0.996224012777970108,
          0.996224012777970108,
          -0.9992829840291237803,
          0.9992829840291237803};
      x = c[n];
    } break;
 
    case 64: {
      T const c[] = {-0.0243502926634244325, 0.0243502926634244325, -0.0729931217877990394, 0.0729931217877990394,
                     -0.1214628192961205544, 0.1214628192961205544, -0.1696444204239928180, 0.1696444204239928180,
                     -0.217423643740007084,  0.217423643740007084,  -0.2646871622087674163, 0.2646871622087674163,
                     -0.3113228719902109561, 0.3113228719902109561, -0.357220158337668115,  0.357220158337668115,
                     -0.4022701579639916036, 0.4022701579639916036, -0.4463660172534640879, 0.4463660172534640879,
                     -0.4894031457070529574, 0.4894031457070529574, -0.531279464019894545,  0.531279464019894545,
                     -0.5718956462026340342, 0.5718956462026340342, -0.6111553551723932502, 0.6111553551723932502,
                     -0.6489654712546573398, 0.6489654712546573398, -0.6852363130542332425, 0.6852363130542332425,
                     -0.7198818501716108268, 0.7198818501716108268, -0.7528199072605318966, 0.7528199072605318966,
                     -0.7839723589433414076, 0.7839723589433414076, -0.8132653151227975597, 0.8132653151227975597,
                     -0.840629296252580362,  0.840629296252580362,  -0.8659993981540928197, 0.8659993981540928197,
                     -0.8893154459951141058, 0.8893154459951141058, -0.9105221370785028057, 0.9105221370785028057,
                     -0.9295691721319395758, 0.9295691721319395758, -0.9464113748584028160, 0.9464113748584028160,
                     -0.9610087996520537189, 0.9610087996520537189, -0.9733268277899109637, 0.9733268277899109637,
                     -0.9833362538846259569, 0.9833362538846259569, -0.9910133714767443207, 0.9910133714767443207,
                     -0.9963401167719552793, 0.9963401167719552793, -0.9993050417357721394, 0.9993050417357721394};
      x           = c[n];
    } break;
  }
  return x;
}
 
template<typename T>
KOKKOS_INLINE_FUNCTION
T
gauss_legendre_weights(Index const m, Index const n)
{
  T x = 0.0;
  switch (m) {
    default: break;
    case 1: {
      T const c[] = {2.000000000000000000};
      x           = c[n];
    } break;
    case 2: {
      T const c[] = {1.000000000000000000, 1.000000000000000000};
      x           = c[n];
    } break;
 
    case 3: {
      T const c[] = {0.8888888888888888888, 0.5555555555555555555, 0.5555555555555555555};
      x           = c[n];
    } break;
 
    case 4: {
      T const c[] = {0.6521451548625461426, 0.6521451548625461426, 0.3478548451374538573, 0.3478548451374538573};
      x           = c[n];
    } break;
 
    case 5: {
      T const c[] = {
          0.5688888888888888888,
          0.4786286704993664680,
          0.4786286704993664680,
          0.2369268850561890875,
          0.2369268850561890875};
      x = c[n];
    } break;
 
    case 6: {
      T const c[] = {
          0.3607615730481386075698335138377,
          0.3607615730481386075698335138377,
          0.4679139345726910473898703439895,
          0.4679139345726910473898703439895,
          0.1713244923791703450402961421727,
          0.1713244923791703450402961421727};
      x = c[n];
    } break;
 
    case 7: {
      T const c[] = {
          0.4179591836734693877,
          0.3818300505051189449503697754889,
          0.3818300505051189449503697754889,
          0.2797053914892766679014677714237,
          0.2797053914892766679014677714237,
          0.1294849661688696932706114326790,
          0.1294849661688696932706114326790};
      x = c[n];
    } break;
 
    case 8: {
      T const c[] = {
          0.3626837833783619829651504492771,
          0.3626837833783619829651504492771,
          0.3137066458778872873379622019866,
          0.3137066458778872873379622019866,
          0.2223810344533744705443559944262,
          0.2223810344533744705443559944262,
          0.1012285362903762591525313543099,
          0.1012285362903762591525313543099};
      x = c[n];
    } break;
 
    case 9: {
      T const c[] = {
          0.3302393550012597631,
          0.1806481606948574040584720312429,
          0.1806481606948574040584720312429,
          0.081274388361574411971892158110523,
          0.081274388361574411971892158110523,
          0.3123470770400028400686304065844,
          0.3123470770400028400686304065844,
          0.2606106964029354623187428694186,
          0.2606106964029354623187428694186};
      x = c[n];
    } break;
 
    case 10: {
      T const c[] = {
          0.295524224714752870,
          0.295524224714752870,
          0.2692667193099963550,
          0.2692667193099963550,
          0.2190863625159820439,
          0.2190863625159820439,
          0.1494513491505805931,
          0.1494513491505805931,
          0.066671344308688137,
          0.066671344308688137};
      x = c[n];
    } break;
 
    case 11: {
      T const c[] = {
          0.2729250867779006307,
          0.2628045445102466621,
          0.2628045445102466621,
          0.2331937645919904799,
          0.2331937645919904799,
          0.1862902109277342514,
          0.1862902109277342514,
          0.1255803694649046246,
          0.1255803694649046246,
          0.0556685671161736664,
          0.0556685671161736664};
      x = c[n];
    } break;
 
    case 12: {
      T const c[] = {
          0.2491470458134027850,
          0.2491470458134027850,
          0.2334925365383548087,
          0.2334925365383548087,
          0.203167426723065921,
          0.203167426723065921,
          0.1600783285433462263,
          0.1600783285433462263,
          0.1069393259953184309,
          0.1069393259953184309,
          0.0471753363865118271,
          0.0471753363865118271};
      x = c[n];
    } break;
 
    case 13: {
      T const c[] = {
          0.2325515532308739101,
          0.2262831802628972384,
          0.2262831802628972384,
          0.2078160475368885023,
          0.2078160475368885023,
          0.1781459807619457382,
          0.1781459807619457382,
          0.138873510219787238,
          0.138873510219787238,
          0.0921214998377284479,
          0.0921214998377284479,
          0.0404840047653158795,
          0.0404840047653158795};
      x = c[n];
    } break;
 
    case 14: {
      T const c[] = {
          0.215263853463157790,
          0.215263853463157790,
          0.2051984637212956039,
          0.2051984637212956039,
          0.1855383974779378137,
          0.1855383974779378137,
          0.1572031671581935345,
          0.1572031671581935345,
          0.1215185706879031846,
          0.1215185706879031846,
          0.0801580871597602098,
          0.0801580871597602098,
          0.0351194603317518630,
          0.0351194603317518630};
      x = c[n];
    } break;
 
    case 15: {
      T const c[] = {
          0.2025782419255612728,
          0.1984314853271115764,
          0.1984314853271115764,
          0.1861610000155622110,
          0.1861610000155622110,
          0.1662692058169939335,
          0.1662692058169939335,
          0.1395706779261543144,
          0.1395706779261543144,
          0.1071592204671719350,
          0.1071592204671719350,
          0.0703660474881081247,
          0.0703660474881081247,
          0.0307532419961172683,
          0.0307532419961172683};
      x = c[n];
    } break;
 
    case 16: {
      T const c[] = {
          0.1894506104550684962,
          0.1894506104550684962,
          0.1826034150449235888,
          0.1826034150449235888,
          0.1691565193950025381,
          0.1691565193950025381,
          0.1495959888165767320,
          0.1495959888165767320,
          0.1246289712555338720,
          0.1246289712555338720,
          0.0951585116824927848,
          0.0951585116824927848,
          0.0622535239386478928,
          0.0622535239386478928,
          0.0271524594117540948,
          0.0271524594117540948};
      x = c[n];
    } break;
 
    case 17: {
      T const c[] = {
          0.1794464703562065254,
          0.1765627053669926463,
          0.1765627053669926463,
          0.1680041021564500445,
          0.1680041021564500445,
          0.154045761076810288,
          0.154045761076810288,
          0.1351363684685254732,
          0.1351363684685254732,
          0.1118838471934039710,
          0.1118838471934039710,
          0.0850361483171791808,
          0.0850361483171791808,
          0.055459529373987201,
          0.055459529373987201,
          0.0241483028685479319,
          0.0241483028685479319};
      x = c[n];
    } break;
 
    case 18: {
      T const c[] = {
          0.1691423829631435918,
          0.1691423829631435918,
          0.1642764837458327229,
          0.1642764837458327229,
          0.1546846751262652449,
          0.1546846751262652449,
          0.1406429146706506512,
          0.1406429146706506512,
          0.1225552067114784601,
          0.1225552067114784601,
          0.1009420441062871655,
          0.1009420441062871655,
          0.0764257302548890565,
          0.0764257302548890565,
          0.0497145488949697964,
          0.0497145488949697964,
          0.0216160135264833103,
          0.0216160135264833103};
      x = c[n];
    } break;
 
    case 19: {
      T const c[] = {
          0.1610544498487836959,
          0.1589688433939543476,
          0.1589688433939543476,
          0.1527660420658596667,
          0.1527660420658596667,
          0.1426067021736066117,
          0.1426067021736066117,
          0.1287539625393362276,
          0.1287539625393362276,
          0.1115666455473339947,
          0.1115666455473339947,
          0.0914900216224499994,
          0.0914900216224499994,
          0.0690445427376412265,
          0.0690445427376412265,
          0.0448142267656996003,
          0.0448142267656996003,
          0.0194617882297264770,
          0.0194617882297264770};
      x = c[n];
    } break;
 
    case 20: {
      T const c[] = {0.1527533871307258506, 0.1527533871307258506, 0.1491729864726037467, 0.1491729864726037467,
                     0.1420961093183820513, 0.1420961093183820513, 0.131688638449176626,  0.131688638449176626,
                     0.1181945319615184173, 0.1181945319615184173, 0.1019301198172404350, 0.1019301198172404350,
                     0.0832767415767047487, 0.0832767415767047487, 0.0626720483341090635, 0.0626720483341090635,
                     0.0406014298003869413, 0.0406014298003869413, 0.017614007139152118,  0.017614007139152118};
      x           = c[n];
    } break;
 
    case 21: {
      T const c[] = {0.1460811336496904271, 0.1445244039899700590, 0.1445244039899700590, 0.1398873947910731547,
                     0.1398873947910731547, 0.1322689386333374617, 0.1322689386333374617, 0.12183141605372853,
                     0.12183141605372853,   0.1087972991671483776, 0.1087972991671483776, 0.0934444234560338615,
                     0.0934444234560338615, 0.0761001136283793020, 0.0761001136283793020, 0.0571344254268572082,
                     0.0571344254268572082, 0.0369537897708524937, 0.0369537897708524937, 0.0160172282577743333,
                     0.0160172282577743333};
      x           = c[n];
    } break;
 
    case 22: {
      T const c[] = {0.139251872855631993,  0.139251872855631993,  0.1365414983460151713, 0.1365414983460151713,
                     0.1311735047870623707, 0.1311735047870623707, 0.1232523768105124242, 0.1232523768105124242,
                     0.1129322960805392183, 0.1129322960805392183, 0.1004141444428809649, 0.1004141444428809649,
                     0.0859416062170677274, 0.0859416062170677274, 0.0697964684245204880, 0.0697964684245204880,
                     0.0522933351526832859, 0.0522933351526832859, 0.0337749015848141547, 0.0337749015848141547,
                     0.0146279952982722006, 0.0146279952982722006};
      x           = c[n];
    } break;
 
    case 23: {
      T const c[] = {0.1336545721861061753, 0.1324620394046966173, 0.1324620394046966173, 0.1289057221880821499,
                     0.1289057221880821499, 0.1230490843067295304, 0.1230490843067295304, 0.1149966402224113649,
                     0.1149966402224113649, 0.1048920914645414100, 0.1048920914645414100, 0.0929157660600351474,
                     0.0929157660600351474, 0.0792814117767189549, 0.0792814117767189549, 0.0642324214085258521,
                     0.0642324214085258521, 0.0480376717310846685, 0.0480376717310846685, 0.0309880058569794443,
                     0.0309880058569794443, 0.0134118594871417720, 0.0134118594871417720};
      x           = c[n];
    } break;
 
    case 24: {
      T const c[] = {0.1279381953467521569, 0.1279381953467521569, 0.1258374563468282961, 0.1258374563468282961,
                     0.1216704729278033912, 0.1216704729278033912, 0.1155056680537256013, 0.1155056680537256013,
                     0.107444270115965634,  0.107444270115965634,  0.097618652104113888,  0.097618652104113888,
                     0.0861901615319532759, 0.0861901615319532759, 0.0733464814110803057, 0.0733464814110803057,
                     0.0592985849154367807, 0.0592985849154367807, 0.0442774388174198061, 0.0442774388174198061,
                     0.0285313886289336631, 0.0285313886289336631, 0.0123412297999871995, 0.0123412297999871995};
      x           = c[n];
    } break;
 
    case 25: {
      T const c[] = {0.1231760537267154512, 0.1222424429903100416, 0.1222424429903100416, 0.1194557635357847722,
                     0.1194557635357847722, 0.1148582591457116483, 0.1148582591457116483, 0.1085196244742636531,
                     0.1085196244742636531, 0.1005359490670506442, 0.1005359490670506442, 0.0910282619829636498,
                     0.0910282619829636498, 0.0801407003350010180, 0.0801407003350010180, 0.0680383338123569172,
                     0.0680383338123569172, 0.0549046959758351919, 0.0549046959758351919, 0.0409391567013063126,
                     0.0409391567013063126, 0.0263549866150321372, 0.0263549866150321372, 0.0113937985010262879,
                     0.0113937985010262879};
      x           = c[n];
    } break;
 
    case 26: {
      T const c[] = {0.1183214152792622765, 0.1183214152792622765, 0.116660443485296582,  0.116660443485296582,
                     0.1133618165463196665, 0.1133618165463196665, 0.1084718405285765906, 0.1084718405285765906,
                     0.102059161094425423,  0.102059161094425423,  0.0942138003559141484, 0.0942138003559141484,
                     0.0850458943134852392, 0.0850458943134852392, 0.0746841497656597458, 0.0746841497656597458,
                     0.0632740463295748355, 0.0632740463295748355, 0.0509758252971478119, 0.0509758252971478119,
                     0.0379623832943627639, 0.0379623832943627639, 0.0244178510926319087, 0.0244178510926319087,
                     0.010551372617343007,  0.010551372617343007};
      x           = c[n];
    } break;
 
    case 27: {
      T const c[] = {0.1142208673789569890, 0.1134763461089651486, 0.1134763461089651486, 0.1112524883568451926,
                     0.1112524883568451926, 0.107578285788533187,  0.107578285788533187,  0.1025016378177457986,
                     0.1025016378177457986, 0.0960887273700285075, 0.0960887273700285075, 0.0884231585437569501,
                     0.0884231585437569501, 0.0796048677730577712, 0.0796048677730577712, 0.0697488237662455929,
                     0.0697488237662455929, 0.058983536859833599,  0.058983536859833599,  0.0474494125206150627,
                     0.0474494125206150627, 0.0352970537574197110, 0.0352970537574197110, 0.0226862315961806231,
                     0.0226862315961806231, 0.0097989960512943602, 0.0097989960512943602};
      x           = c[n];
    } break;
 
    case 28: {
      T const c[] = {0.1100470130164751962, 0.1100470130164751962, 0.108711192258294135,  0.108711192258294135,
                     0.1060557659228464179, 0.1060557659228464179, 0.1021129675780607698, 0.1021129675780607698,
                     0.0969306579979299158, 0.0969306579979299158, 0.0905717443930328409, 0.0905717443930328409,
                     0.08311341722890121,   0.08311341722890121,   0.0746462142345687790, 0.0746462142345687790,
                     0.0652729239669995957, 0.0652729239669995957, 0.0551073456757167454, 0.0551073456757167454,
                     0.0442729347590042278, 0.0442729347590042278, 0.0329014277823043799, 0.0329014277823043799,
                     0.0211321125927712597, 0.0211321125927712597, 0.009124282593094517,  0.009124282593094517};
      x           = c[n];
    } break;
 
    case 29: {
      T const c[] = {0.1064793817183142442, 0.1058761550973209414, 0.1058761550973209414, 0.1040733100777293739,
                     0.1040733100777293739, 0.1010912737599149661, 0.1010912737599149661, 0.0969638340944086063,
                     0.0969638340944086063, 0.091737757139258763,  0.091737757139258763,  0.0854722573661725275,
                     0.0854722573661725275, 0.0782383271357637838, 0.0782383271357637838, 0.0701179332550512785,
                     0.0701179332550512785, 0.0612030906570791385, 0.0612030906570791385, 0.051594826902497923,
                     0.051594826902497923,  0.0414020625186828361, 0.0414020625186828361, 0.0307404922020936226,
                     0.0307404922020936226, 0.0197320850561227059, 0.0197320850561227059, 0.0085169038787464096,
                     0.0085169038787464096};
      x           = c[n];
    } break;
 
    case 30: {
      T const c[] = {0.1028526528935588403, 0.1028526528935588403, 0.1017623897484055045, 0.1017623897484055045,
                     0.09959342058679526,   0.09959342058679526,   0.0963687371746442596, 0.0963687371746442596,
                     0.0921225222377861287, 0.0921225222377861287, 0.0868997872010829798, 0.0868997872010829798,
                     0.0807558952294202153, 0.0807558952294202153, 0.073755974737705206,  0.073755974737705206,
                     0.0659742298821804951, 0.0659742298821804951, 0.0574931562176190664, 0.0574931562176190664,
                     0.0484026728305940529, 0.0484026728305940529, 0.0387991925696270495, 0.0387991925696270495,
                     0.0287847078833233693, 0.0287847078833233693, 0.0184664683110909591, 0.0184664683110909591,
                     0.0079681924961666056, 0.0079681924961666056};
      x           = c[n];
    } break;
 
    case 31: {
      T const c[] = {0.099720544793426451,  0.0992250112266723078, 0.0992250112266723078, 0.0977433353863287250,
                     0.0977433353863287250, 0.0952902429123195128, 0.0952902429123195128, 0.0918901138936414782,
                     0.0918901138936414782, 0.0875767406084778761, 0.0875767406084778761, 0.0823929917615892639,
                     0.0823929917615892639, 0.0763903865987766164, 0.0763903865987766164, 0.0696285832354103661,
                     0.0696285832354103661, 0.0621747865610284269, 0.0621747865610284269, 0.0541030824249168537,
                     0.0541030824249168537, 0.0454937075272011029, 0.0454937075272011029, 0.0364322739123854640,
                     0.0364322739123854640, 0.0270090191849794218, 0.0270090191849794218, 0.0173186207903105824,
                     0.0173186207903105824, 0.007470831579248775,  0.007470831579248775};
      x           = c[n];
    } break;
 
    case 32: {
      T const c[] = {0.0965400885147278005, 0.0965400885147278005, 0.0956387200792748594, 0.0956387200792748594,
                     0.0938443990808045656, 0.0938443990808045656, 0.0911738786957638847, 0.0911738786957638847,
                     0.0876520930044038111, 0.0876520930044038111, 0.0833119242269467552, 0.0833119242269467552,
                     0.078193895787070306,  0.078193895787070306,  0.0723457941088485062, 0.0723457941088485062,
                     0.0658222227763618468, 0.0658222227763618468, 0.0586840934785355471, 0.0586840934785355471,
                     0.0509980592623761761, 0.0509980592623761761, 0.0428358980222266806, 0.0428358980222266806,
                     0.0342738629130214331, 0.0342738629130214331, 0.0253920653092620594, 0.0253920653092620594,
                     0.0162743947309056706, 0.0162743947309056706, 0.0070186100094700966, 0.0070186100094700966};
      x           = c[n];
    } break;
 
    case 33: {
      T const c[] = {0.0937684461602099965, 0.093356426065596116,  0.093356426065596116,  0.0921239866433168462,
                     0.0921239866433168462, 0.0900819586606385772, 0.0900819586606385772, 0.0872482876188443376,
                     0.0872482876188443376, 0.0836478760670387076, 0.0836478760670387076, 0.0793123647948867383,
                     0.0793123647948867383, 0.0742798548439541493, 0.0742798548439541493, 0.0685945728186567128,
                     0.0685945728186567128, 0.0623064825303174800, 0.0623064825303174800, 0.0554708466316635612,
                     0.0554708466316635612, 0.0481477428187116956, 0.0481477428187116956, 0.0404015413316695915,
                     0.0404015413316695915, 0.0323003586323289532, 0.0323003586323289532, 0.0239155481017494803,
                     0.0239155481017494803, 0.0153217015129346761, 0.0153217015129346761, 0.0066062278475873780,
                     0.0066062278475873780};
      x           = c[n];
    } break;
 
    case 34: {
      T const c[] = {0.0909567403302598736, 0.0909567403302598736, 0.0902030443706407295, 0.0902030443706407295,
                     0.0887018978356938692, 0.0887018978356938692, 0.0864657397470357497, 0.0864657397470357497,
                     0.0835130996998456551, 0.0835130996998456551, 0.0798684443397718447, 0.0798684443397718447,
                     0.0755619746600319312, 0.0755619746600319312, 0.0706293758142557249, 0.0706293758142557249,
                     0.0651115215540764113, 0.0651115215540764113, 0.059054135827524493,  0.059054135827524493,
                     0.0525074145726781061, 0.0525074145726781061, 0.045525611523353272,  0.045525611523353272,
                     0.038166593796387516,  0.038166593796387516,  0.0304913806384461318, 0.0304913806384461318,
                     0.0225637219854949700, 0.0225637219854949700, 0.014450162748595035,  0.014450162748595035,
                     0.0062291405559086847, 0.0062291405559086847};
      x           = c[n];
    } break;
 
    case 35: {
      T const c[] = {0.0884867949071042906, 0.0881405304302754629, 0.0881405304302754629, 0.0871044469971835342,
                     0.0871044469971835342, 0.0853866533920991252, 0.0853866533920991252, 0.0830005937288565883,
                     0.0830005937288565883, 0.0799649422423242629, 0.0799649422423242629, 0.0763034571554420535,
                     0.0763034571554420535, 0.0720447947725600646, 0.0720447947725600646, 0.0672222852690869039,
                     0.0672222852690869039, 0.0618736719660801888, 0.0618736719660801888, 0.0560408162123701285,
                     0.0560408162123701285, 0.0497693704013535298, 0.0497693704013535298, 0.0431084223261702187,
                     0.0431084223261702187, 0.0361101158634633805, 0.0361101158634633805, 0.0288292601088942540,
                     0.0288292601088942540, 0.0213229799114835808, 0.0213229799114835808, 0.0136508283483614922,
                     0.0136508283483614922, 0.0058834334204430849, 0.0058834334204430849};
      x           = c[n];
    } break;
 
    case 36: {
      T const c[] = {0.0859832756703947474, 0.0859832756703947474, 0.0853466857393386274, 0.0853466857393386274,
                     0.0840782189796619349, 0.0840782189796619349, 0.0821872667043397095, 0.0821872667043397095,
                     0.0796878289120716019, 0.0796878289120716019, 0.07659841064587067,   0.07659841064587067,
                     0.0729418850056530613, 0.0729418850056530613, 0.0687453238357364426, 0.0687453238357364426,
                     0.0640397973550154895, 0.0640397973550154895, 0.0588601442453248173, 0.0588601442453248173,
                     0.0532447139777599190, 0.0532447139777599190, 0.047235083490265978,  0.047235083490265978,
                     0.0408757509236448954, 0.0408757509236448954, 0.0342138107703072299, 0.0342138107703072299,
                     0.027298621498568779,  0.027298621498568779,  0.0201815152977354715, 0.0201815152977354715,
                     0.0129159472840655744, 0.0129159472840655744, 0.0055657196642450453, 0.0055657196642450453};
      x           = c[n];
    } break;
 
    case 37: {
      T const c[] = {0.0837683609931389047, 0.083474573625862787,  0.083474573625862787,  0.082595272236437250,
                     0.082595272236437250,  0.0811366245084650305, 0.0811366245084650305, 0.0791088618375293807,
                     0.0791088618375293807, 0.0765262075705292378, 0.0765262075705292378, 0.073406777248488172,
                     0.073406777248488172,  0.0697724515557003448, 0.0697724515557003448, 0.0656487228727512494,
                     0.0656487228727512494, 0.061064516523225986,  0.061064516523225986,  0.0560519879982749178,
                     0.0560519879982749178, 0.05064629765482460,   0.05064629765482460,   0.0448853646624371666,
                     0.0448853646624371666, 0.0388096025019345444, 0.0388096025019345444, 0.0324616398475214810,
                     0.0324616398475214810, 0.0258860369905589335, 0.0258860369905589335, 0.0191290444890839660,
                     0.0191290444890839660, 0.0122387801003075565, 0.0122387801003075565, 0.0052730572794979393,
                     0.0052730572794979393};
      x           = c[n];
    } break;
 
    case 38: {
      T const c[] = {0.0815250292803857866, 0.0815250292803857866, 0.0809824937705971006, 0.0809824937705971006,
                     0.0799010332435278215, 0.0799010332435278215, 0.0782878446582109480, 0.0782878446582109480,
                     0.0761536635484463960, 0.0761536635484463960, 0.0735126925847434571, 0.0735126925847434571,
                     0.0703825070668989547, 0.0703825070668989547, 0.0667839379791404119, 0.0667839379791404119,
                     0.062740933392133054,  0.062740933392133054,  0.0582803991469972060, 0.0582803991469972060,
                     0.0534320199103323199, 0.0534320199103323199, 0.0482280618607586833, 0.0482280618607586833,
                     0.0427031585046744342, 0.0427031585046744342, 0.0368940815940247381, 0.0368940815940247381,
                     0.0308395005451750546, 0.0308395005451750546, 0.0245797397382323758, 0.0245797397382323758,
                     0.0181565777096132368, 0.0181565777096132368, 0.0116134447164686741, 0.0116134447164686741,
                     0.0050028807496393456, 0.0050028807496393456};
      x           = c[n];
    } break;
 
    case 39: {
      T const c[] = {0.0795276221394428524, 0.0792762225683684710, 0.0792762225683684710, 0.0785236132873711767,
                     0.0785236132873711767, 0.0772745525446820167, 0.0772745525446820167, 0.075536937322836057,
                     0.075536937322836057,  0.0733217534142686173, 0.0733217534142686173, 0.0706430059706087607,
                     0.0706430059706087607, 0.0675176309662312653, 0.0675176309662312653, 0.0639653881386823889,
                     0.0639653881386823889, 0.0600087360885961495, 0.0600087360885961495, 0.0556726903409162999,
                     0.0556726903409162999, 0.05098466529212940,   0.05098466529212940,   0.0459743011089166318,
                     0.0459743011089166318, 0.040673276847933843,  0.040673276847933843,  0.0351151114981313307,
                     0.0351151114981313307, 0.0293349559839033785, 0.0293349559839033785, 0.023369384832178164,
                     0.023369384832178164,  0.0172562290937249190, 0.0172562290937249190, 0.0110347889391645942,
                     0.0110347889391645942, 0.0047529446916351013, 0.0047529446916351013};
      x           = c[n];
    } break;
 
    case 40: {
      T const c[] = {0.0775059479784248112, 0.0775059479784248112, 0.0770398181642479655, 0.0770398181642479655,
                     0.0761103619006262423, 0.0761103619006262423, 0.0747231690579682642, 0.0747231690579682642,
                     0.0728865823958040590, 0.0728865823958040590, 0.0706116473912867796, 0.0706116473912867796,
                     0.0679120458152339038, 0.0679120458152339038, 0.0648040134566010380, 0.0648040134566010380,
                     0.0613062424929289391, 0.0613062424929289391, 0.0574397690993915513, 0.0574397690993915513,
                     0.0532278469839368243, 0.0532278469839368243, 0.0486958076350722320, 0.0486958076350722320,
                     0.0438709081856732719, 0.0438709081856732719, 0.03878216797447201,   0.03878216797447201,
                     0.0334601952825478473, 0.0334601952825478473, 0.0279370069800234010, 0.0279370069800234010,
                     0.0222458491941669572, 0.0222458491941669572, 0.0164210583819078887, 0.0164210583819078887,
                     0.0104982845311528136, 0.0104982845311528136, 0.004521277098533191,  0.004521277098533191};
      x           = c[n];
    } break;
 
    case 41: {
      T const c[] = {0.0756955356472983723, 0.0754787470927158240, 0.0754787470927158240, 0.0748296231762215518,
                     0.0748296231762215518, 0.0737518820272234699, 0.0737518820272234699, 0.0722516968610230733,
                     0.0722516968610230733, 0.0703376606208174974, 0.0703376606208174974, 0.0680207367608767667,
                     0.0680207367608767667, 0.065314196453527410,  0.065314196453527410,  0.0622335425809663164,
                     0.0622335425809663164, 0.0587964209498719449, 0.0587964209498719449, 0.0550225192425787418,
                     0.0550225192425787418, 0.0509334542946174947, 0.0509334542946174947, 0.0465526483690143420,
                     0.0465526483690143420, 0.0419051951959096894, 0.0419051951959096894, 0.0370177167035079884,
                     0.0370177167035079884, 0.0319182117316992817, 0.0319182117316992817, 0.0266358992071104454,
                     0.0266358992071104454, 0.0212010633687795530, 0.0212010633687795530, 0.0156449384078185885,
                     0.0156449384078185885, 0.009999938773905945,  0.009999938773905945,  0.0043061403581648876,
                     0.0043061403581648876};
      x           = c[n];
    } break;
 
    case 42: {
      T const c[] = {0.0738642342321728799, 0.0738642342321728799, 0.073460813453467528,  0.073460813453467528,
                     0.0726561752438041048, 0.0726561752438041048, 0.0714547142651709829, 0.0714547142651709829,
                     0.0698629924925941597, 0.0698629924925941597, 0.0678897033765219448, 0.0678897033765219448,
                     0.0655456243649089789, 0.0655456243649089789, 0.0628435580450025764, 0.0628435580450025764,
                     0.0597982622275866543, 0.0597982622275866543, 0.0564263693580183816, 0.0564263693580183816,
                     0.0527462956991740703, 0.0527462956991740703, 0.0487781407928032450, 0.0487781407928032450,
                     0.04454357777196587,   0.04454357777196587,   0.0400657351806922617, 0.0400657351806922617,
                     0.03536907109759211,   0.03536907109759211,   0.0304792406996034683, 0.0304792406996034683,
                     0.0254229595261130478, 0.0254229595261130478, 0.0202278695690526447, 0.0202278695690526447,
                     0.01492244369735749,   0.01492244369735749,   0.0095362203017485024, 0.0095362203017485024,
                     0.0041059986046490846, 0.0041059986046490846};
      x           = c[n];
    } break;
 
    case 43: {
      T const c[] = {0.0722157516937989879, 0.0720275019714219743, 0.0720275019714219743, 0.0714637342525141412,
                     0.0714637342525141412, 0.0705273877650850281, 0.0705273877650850281, 0.0692233441936566842,
                     0.0692233441936566842, 0.0675584022293651691, 0.0675584022293651691, 0.0655412421263227974,
                     0.0655412421263227974, 0.0631823804493961123, 0.0631823804493961123, 0.0604941152499912945,
                     0.0604941152499912945, 0.0574904619569105194, 0.0574904619569105194, 0.0541870803188817868,
                     0.0541870803188817868, 0.0506011927843901565, 0.0506011927843901565, 0.0467514947543465800,
                     0.0467514947543465800, 0.0426580571979820837, 0.0426580571979820837, 0.0383422221941326575,
                     0.0383422221941326575, 0.0338264920868602923, 0.0338264920868602923, 0.029134413261498494,
                     0.029134413261498494,  0.0242904566138388159, 0.0242904566138388159, 0.0193199014236839003,
                     0.0193199014236839003, 0.0142487564315764861, 0.0142487564315764861, 0.0091039966374014033,
                     0.0091039966374014033, 0.0039194902538441272, 0.0039194902538441272};
      x           = c[n];
    } break;
 
    case 44: {
      T const c[] = {0.0705491577893540688, 0.0705491577893540688, 0.0701976854735582125, 0.0701976854735582125,
                     0.0694964918615725780, 0.0694964918615725780, 0.0684490702693666609, 0.0684490702693666609,
                     0.0670606389062936523, 0.0670606389062936523, 0.0653381148791814349, 0.0653381148791814349,
                     0.0632900797332038549, 0.0632900797332038549, 0.0609267367015619680, 0.0609267367015619680,
                     0.0582598598775954953, 0.0582598598775954953, 0.0553027355637280525, 0.0553027355637280525,
                     0.0520700960917044618, 0.0520700960917044618, 0.0485780464483520375, 0.0485780464483520375,
                     0.0448439840819700314, 0.0448439840819700314, 0.0408865123103462189, 0.0408865123103462189,
                     0.0367253478138088736, 0.0367253478138088736, 0.0323812228120698208, 0.0323812228120698208,
                     0.0278757828212810100, 0.0278757828212810100, 0.0232314819020192106, 0.0232314819020192106,
                     0.0184714817368147491, 0.0184714817368147491, 0.0136195867555799855, 0.0136195867555799855,
                     0.008700481367524844,  0.008700481367524844,  0.0037454048031127775, 0.0037454048031127775};
      x           = c[n];
    } break;
 
    case 45: {
      T const c[] = {0.0690418248292320201, 0.0688773169776613228, 0.0688773169776613228, 0.0683845773786696745,
                     0.0683845773786696745, 0.067565954163607536,  0.067565954163607536,  0.0664253484498425280,
                     0.0664253484498425280, 0.0649681957507234308, 0.0649681957507234308, 0.0632014400738199377,
                     0.0632014400738199377, 0.0611335008310665225, 0.0611335008310665225, 0.0587742327188417385,
                     0.0587742327188417385, 0.056134878759786476,  0.056134878759786476,  0.0532280167312689519,
                     0.0532280167312689519, 0.0500674992379520297, 0.0500674992379520297, 0.0466683877183733652,
                     0.0466683877183733652, 0.043046880709164971,  0.043046880709164971,  0.0392202367293024475,
                     0.0392202367293024475, 0.0352066922016090162, 0.0352066922016090162, 0.0310253749345154671,
                     0.0310253749345154671, 0.0266962139675776648, 0.0266962139675776648, 0.0222398475505787323,
                     0.0222398475505787323, 0.0176775352579375906, 0.0176775352579375906, 0.0130311049915827843,
                     0.0130311049915827843, 0.0083231892962182416, 0.0083231892962182416, 0.003582663155283558,
                     0.003582663155283558};
      x           = c[n];
    } break;
 
    case 46: {
      T const c[] = {0.0675186858490364588, 0.0675186858490364588, 0.0672106136006781758, 0.0672106136006781758,
                     0.0665958747684548873, 0.0665958747684548873, 0.0656772742677812073, 0.0656772742677812073,
                     0.0644590034671390695, 0.0644590034671390695, 0.0629466210643945081, 0.0629466210643945081,
                     0.0611470277246504810, 0.0611470277246504810, 0.0590684345955463148, 0.0590684345955463148,
                     0.0567203258439912358, 0.0567203258439912358, 0.054113415385856754,  0.054113415385856754,
                     0.0512595980071430213, 0.0512595980071430213, 0.048171895101712200,  0.048171895101712200,
                     0.044864395277318126,  0.044864395277318126,  0.0413521901096787297, 0.0413521901096787297,
                     0.0376513053573860713, 0.0376513053573860713, 0.033778627999106896,  0.033778627999106896,
                     0.0297518295522027557, 0.0297518295522027557, 0.0255892863971300106, 0.0255892863971300106,
                     0.0213099987541365010, 0.0213099987541365010, 0.0169335140078362380, 0.0169335140078362380,
                     0.0124798837709886842, 0.0124798837709886842, 0.0079698982297246224, 0.0079698982297246224,
                     0.0034303008681070482, 0.0034303008681070482};
      x           = c[n];
    } break;
 
    case 47: {
      T const c[] = {0.0661351296236554796, 0.0659905335888104745, 0.0659905335888104745, 0.0655573777665497402,
                     0.0655573777665497402, 0.0648375562389457267, 0.0648375562389457267, 0.0638342166057170306,
                     0.0638342166057170306, 0.0625517462209216626, 0.0625517462209216626, 0.0609957530087396453,
                     0.0609957530087396453, 0.0591730409423388759, 0.0591730409423388759, 0.0570915802932315402,
                     0.0570915802932315402, 0.0547604727815302259, 0.0547604727815302259, 0.0521899117800571448,
                     0.0521899117800571448, 0.0493911377473611696, 0.0493911377473611696, 0.0463763890865059112,
                     0.0463763890865059112, 0.0431588486484795382, 0.0431588486484795382, 0.0397525861225310037,
                     0.0397525861225310037, 0.0361724965841749516, 0.0361724965841749516, 0.0324342355151847567,
                     0.0324342355151847567, 0.0285541507006433865, 0.0285541507006433865, 0.0245492116596588185,
                     0.0245492116596588185, 0.0204369381476684276, 0.0204369381476684276, 0.0162353331464330596,
                     0.0162353331464330596, 0.0119628484643123209, 0.0119628484643123209, 0.0076386162958488336,
                     0.0076386162958488336, 0.0032874538425280148, 0.0032874538425280148};
      x           = c[n];
    } break;
 
    case 48: {
      T const c[] = {0.0647376968126839225, 0.0647376968126839225, 0.0644661644359500822, 0.0644661644359500822,
                     0.0639242385846481866, 0.0639242385846481866, 0.063114192286254025,  0.063114192286254025,
                     0.0620394231598926639, 0.0620394231598926639, 0.060704439165893880,  0.060704439165893880,
                     0.0591148396983956357, 0.0591148396983956357, 0.0572772921004032157, 0.0572772921004032157,
                     0.0551995036999841628, 0.0551995036999841628, 0.0528901894851936670, 0.0528901894851936670,
                     0.0503590355538544749, 0.0503590355538544749, 0.0476166584924904748, 0.0476166584924904748,
                     0.04467456085669428,   0.04467456085669428,   0.0415450829434647492, 0.0415450829434647492,
                     0.0382413510658307063, 0.0382413510658307063, 0.0347772225647704388, 0.0347772225647704388,
                     0.0311672278327980889, 0.0311672278327980889, 0.0274265097083569482, 0.0274265097083569482,
                     0.0235707608393243791, 0.0235707608393243791, 0.0196161604573555278, 0.0196161604573555278,
                     0.0155793157229438487, 0.0155793157229438487, 0.0114772345792345394, 0.0114772345792345394,
                     0.0073275539012762621, 0.0073275539012762621, 0.0031533460523058386, 0.0031533460523058386};
      x           = c[n];
    } break;
 
    case 49: {
      T const c[] = {0.0634632814047905977, 0.0633355092964917485, 0.0633355092964917485, 0.0629527074651956994,
                     0.0629527074651956994, 0.0623164173200572674, 0.0623164173200572674, 0.0614292009791929362,
                     0.0614292009791929362, 0.0602946309531520173, 0.0602946309531520173, 0.058917275760027266,
                     0.058917275760027266,  0.0573026815301874754, 0.0573026815301874754, 0.0554573496748035886,
                     0.0554573496748035886, 0.0533887107082589685, 0.0533887107082589685, 0.0511050943301445906,
                     0.0511050943301445906, 0.0486156958878282402, 0.0486156958878282402, 0.0459305393555958535,
                     0.0459305393555958535, 0.0430604369812595979, 0.0430604369812595979, 0.0400169457663730213,
                     0.0400169457663730213, 0.0368123209630006898, 0.0368123209630006898, 0.0334594667916221743,
                     0.0334594667916221743, 0.0299718846205838253, 0.0299718846205838253, 0.0263636189270660169,
                     0.0263636189270660169, 0.0226492015874466764, 0.0226492015874466764, 0.0188435958530894584,
                     0.0188435958530894584, 0.014962144935624651,  0.014962144935624651,  0.0110205510315935804,
                     0.0110205510315935804, 0.0070350995900864514, 0.0070350995900864514, 0.0030272789889229050,
                     0.0030272789889229050};
      x           = c[n];
    } break;
 
    case 50: {
      T const c[] = {0.0621766166553472623, 0.0621766166553472623, 0.0619360674206832433, 0.0619360674206832433,
                     0.0614558995903166637, 0.0614558995903166637, 0.0607379708417702160, 0.0607379708417702160,
                     0.0597850587042654575, 0.0597850587042654575, 0.0586008498132224458, 0.0586008498132224458,
                     0.0571899256477283837, 0.0571899256477283837, 0.0555577448062125176, 0.0555577448062125176,
                     0.0537106218889962465, 0.0537106218889962465, 0.051655703069581138,  0.051655703069581138,
                     0.0494009384494663149, 0.0494009384494663149, 0.04695505130394843,   0.04695505130394843,
                     0.0443275043388032754, 0.0443275043388032754, 0.0415284630901476974, 0.0415284630901476974,
                     0.038568756612587675,  0.038568756612587675,  0.0354598356151461541, 0.0354598356151461541,
                     0.0322137282235780166, 0.0322137282235780166, 0.0288429935805351980, 0.0288429935805351980,
                     0.0253606735700123904, 0.0253606735700123904, 0.021780243170124792,  0.021780243170124792,
                     0.0181155607134893903, 0.0181155607134893903, 0.0143808227614855744, 0.0143808227614855744,
                     0.0105905483836509692, 0.0105905483836509692, 0.0067597991957454015, 0.0067597991957454015,
                     0.0029086225531551409, 0.0029086225531551409};
      x           = c[n];
    } break;
 
    case 51: {
      T const c[] = {0.0609989248412058801, 0.0608854648448563438, 0.0608854648448563438, 0.0605455069347377951,
                     0.0605455069347377951, 0.0599803157775032520, 0.0599803157775032520, 0.0591919939229615437,
                     0.0591919939229615437, 0.0581834739825921405, 0.0581834739825921405, 0.0569585077202586621,
                     0.0569585077202586621, 0.0555216520957386930, 0.0555216520957386930, 0.0538782523130455614,
                     0.0538782523130455614, 0.0520344219366970875, 0.0520344219366970875, 0.0499970201500574097,
                     0.0499970201500574097, 0.0477736262406231019, 0.0477736262406231019, 0.0453725114076500687,
                     0.0453725114076500687, 0.0428026079978800866, 0.0428026079978800866, 0.0400734762854964531,
                     0.0400734762854964531, 0.0371952689232602928, 0.0371952689232602928, 0.0341786932041883362,
                     0.0341786932041883362, 0.0310349712901600084, 0.0310349712901600084, 0.0277757985941624771,
                     0.0277757985941624771, 0.0244133005737814342, 0.0244133005737814342, 0.0209599884017032105,
                     0.0209599884017032105, 0.0174287147234010522, 0.0174287147234010522, 0.0138326340064778222,
                     0.0138326340064778222, 0.0101851912978217299, 0.0101851912978217299, 0.0065003377832526002,
                     0.0065003377832526002, 0.0027968071710898955, 0.0027968071710898955};
      x           = c[n];
    } break;
 
    case 52: {
      T const c[] = {0.0598103657452918602, 0.0598103657452918602, 0.0595962601712481582, 0.0595962601712481582,
                     0.0591688154660429703, 0.0591688154660429703, 0.0585295617718138685, 0.0585295617718138685,
                     0.0576807874525268276, 0.0576807874525268276, 0.0566255309023685971, 0.0566255309023685971,
                     0.0553675696693026525, 0.0553675696693026525, 0.0539114069327572647, 0.0539114069327572647,
                     0.0522622553839069930, 0.0522622553839069930, 0.0504260185663423772, 0.0504260185663423772,
                     0.0484092697440748968, 0.0484092697440748968, 0.046219228372784793,  0.046219228372784793,
                     0.0438637342590004079, 0.0438637342590004079, 0.0413512195005602716, 0.0413512195005602716,
                     0.0386906783104239789, 0.0386906783104239789, 0.0358916348350972329, 0.0358916348350972329,
                     0.0329641090897187979, 0.0329641090897187979, 0.0299185811471439466, 0.0299185811471439466,
                     0.0267659537465040134, 0.0267659537465040134, 0.0235175135539844615, 0.0235175135539844615,
                     0.0201848915079807922, 0.0201848915079807922, 0.0167800233963007356, 0.0167800233963007356,
                     0.0133151149823409606, 0.0133151149823409606, 0.0098026345794627520, 0.0098026345794627520,
                     0.0062555239629732768, 0.0062555239629732768, 0.0026913169500471111, 0.0026913169500471111};
      x           = c[n];
    } break;
 
    case 53: {
      T const c[] = {0.0587187941511643645, 0.0586175862327202633, 0.0586175862327202633, 0.0583143113622560075,
                     0.0583143113622560075, 0.0578100149917131963, 0.0578100149917131963, 0.0571064355362671917,
                     0.0571064355362671917, 0.0562059983817397098, 0.0562059983817397098, 0.0551118075239335990,
                     0.0551118075239335990, 0.0538276348687310290, 0.0538276348687310290, 0.052357907229872718,
                     0.052357907229872718,  0.0507076910692927152, 0.0507076910692927152, 0.0488826750326991404,
                     0.0488826750326991404, 0.0468891503407503140, 0.0468891503407503140, 0.0447339891036728102,
                     0.0447339891036728102, 0.0424246206345200135, 0.0424246206345200135, 0.0399690058435403821,
                     0.0399690058435403821, 0.0373756098034829156, 0.0373756098034829156, 0.0346533725835342379,
                     0.0346533725835342379, 0.0318116784590193230, 0.0318116784590193230, 0.0288603236178237362,
                     0.0288603236178237362, 0.0258094825107575177, 0.0258094825107575177, 0.0226696730570702083,
                     0.0226696730570702083, 0.0194517211076368953, 0.0194517211076368953, 0.016166725256687463,
                     0.016166725256687463,  0.0128260261442403791, 0.0128260261442403791, 0.0094412022849403443,
                     0.0094412022849403443, 0.0060242762269486732, 0.0060242762269486732, 0.0025916837205670318,
                     0.0025916837205670318};
      x           = c[n];
    } break;
 
    case 54: {
      T const c[] = {0.0576175367071470246, 0.0576175367071470246, 0.0574261370541121148, 0.0574261370541121148,
                     0.0570439735587945985, 0.0570439735587945985, 0.0564723157306259650, 0.0564723157306259650,
                     0.0557130625605899876, 0.0557130625605899876, 0.0547687362130579863, 0.0547687362130579863,
                     0.0536424736475536112, 0.0536424736475536112, 0.0523380161982987446, 0.0523380161982987446,
                     0.0508596971461881443, 0.0508596971461881443, 0.0492124273245288860, 0.0492124273245288860,
                     0.0474016788064449910, 0.0474016788064449910, 0.0454334667282767139, 0.0454334667282767139,
                     0.0433143293095970154, 0.0433143293095970154, 0.0410513061366449742, 0.0410513061366449742,
                     0.0386519147821025168, 0.0386519147821025168, 0.0361241258403835525, 0.0361241258403835525,
                     0.0334763364643726457, 0.0334763364643726457, 0.030717342497870676,  0.030717342497870676,
                     0.0278563093105958702, 0.0278563093105958702, 0.0249027414672087730, 0.0249027414672087730,
                     0.0218664514228530859, 0.0218664514228530859, 0.0187575276214693779, 0.0187575276214693779,
                     0.015586303035924131,  0.015586303035924131,  0.0123633281288476441, 0.0123633281288476441,
                     0.0090993694555093969, 0.0090993694555093969, 0.0058056110152399848, 0.0058056110152399848,
                     0.0024974818357615857, 0.0024974818357615857};
      x           = c[n];
    } break;
 
    case 55: {
      T const c[] = {0.0566029764445604254, 0.056512318249772001,  0.056512318249772001,  0.0562406340710843680,
                     0.0562406340710843680, 0.0557887941952840871, 0.0557887941952840871, 0.055158246002508687,
                     0.055158246002508687,  0.054351009329911102,  0.054351009329911102,  0.0533696700016054727,
                     0.0533696700016054727, 0.052217371545632084,  0.052217371545632084,  0.0508978051244939792,
                     0.0508978051244939792, 0.049415197711551739,  0.049415197711551739,  0.0477742985512006955,
                     0.0477742985512006955, 0.0459803639462838381, 0.0459803639462838381, 0.0440391404216065898,
                     0.0440391404216065898, 0.041956846317718762,  0.041956846317718762,  0.0397401518743371796,
                     0.0397401518743371796, 0.0373961578679655452, 0.0373961578679655452, 0.0349323728735898874,
                     0.0349323728735898874, 0.0323566892261858316, 0.0323566892261858316, 0.0296773577651610412,
                     0.0296773577651610412, 0.0269029614563962706, 0.0269029614563962706, 0.0240423880097256220,
                     0.0240423880097256220, 0.0211048016680164541, 0.0211048016680164541, 0.0180996145207290624,
                     0.0180996145207290624, 0.0150364583335117882, 0.0150364583335117882, 0.0119251607198486121,
                     0.0119251607198486121, 0.0087757461070585281, 0.0087757461070585281, 0.0055986322665607673,
                     0.0055986322665607673, 0.0024083236199797888, 0.0024083236199797888};
      x           = c[n];
    } break;
 
    case 56: {
      T const c[] = {0.0555797463065143958, 0.0555797463065143958, 0.0554079525032451232, 0.0554079525032451232,
                     0.0550648959017624257, 0.0550648959017624257, 0.0545516368708894210, 0.0545516368708894210,
                     0.0538697618657144857, 0.0538697618657144857, 0.0530213785240107639, 0.0530213785240107639,
                     0.0520091091517413998, 0.0520091091517413998, 0.0508360826177984805, 0.0508360826177984805,
                     0.0495059246830475789, 0.0495059246830475789, 0.0480227467936002581, 0.0480227467936002581,
                     0.0463911333730018967, 0.0463911333730018967, 0.0446161276526922832, 0.0446161276526922832,
                     0.0427032160846670865, 0.0427032160846670865, 0.0406583113847445178, 0.0406583113847445178,
                     0.0384877342592476624, 0.0384877342592476624, 0.0361981938723151860, 0.0361981938723151860,
                     0.0337967671156117612, 0.0337967671156117612, 0.0312908767473104478, 0.0312908767473104478,
                     0.0286882684738227417, 0.0286882684738227417, 0.0259969870583919521, 0.0259969870583919521,
                     0.023225351562565316,  0.023225351562565316,  0.0203819298824025726, 0.0203819298824025726,
                     0.0174755129114009465, 0.0174755129114009465, 0.0145150892780214718, 0.0145150892780214718,
                     0.0115098243403833821, 0.0115098243403833821, 0.0084690631633078876, 0.0084690631633078876,
                     0.005402522246015337,  0.005402522246015337,  0.0023238553757732155, 0.0023238553757732155};
      x           = c[n];
    } break;
 
    case 57: {
      T const c[] = {0.0546343287565840240, 0.0545528036047618864, 0.0545528036047618864, 0.0543084714524986431,
                     0.0543084714524986431, 0.0539020614832985746, 0.0539020614832985746, 0.0533347865848191584,
                     0.0533347865848191584, 0.052608339729177432,  0.052608339729177432,  0.0517248889205178247,
                     0.0517248889205178247, 0.0506870707249274086, 0.0506870707249274086, 0.0494979824020196789,
                     0.0494979824020196789, 0.0481611726616877512, 0.0481611726616877512, 0.0466806310736415037,
                     0.0466806310736415037, 0.0450607761613811577, 0.0450607761613811577, 0.043306442216215196,
                     0.043306442216215196,  0.0414228648708011103, 0.0414228648708011103, 0.0394156654754801140,
                     0.0394156654754801140, 0.0372908343244173173, 0.0372908343244173173, 0.0350547127823126175,
                     0.0350547127823126175, 0.032713974366371568,  0.032713974366371568,  0.0302756048426939994,
                     0.0302756048426939994, 0.0277468814021801923, 0.0277468814021801923, 0.0251353509909181226,
                     0.0251353509909181226, 0.0224488078907764380, 0.0224488078907764380, 0.0196952706994885203,
                     0.0196952706994885203, 0.0168829590234415490, 0.0168829590234415490, 0.0140202707907535561,
                     0.0140202707907535561, 0.0111157637323359901, 0.0111157637323359901, 0.0081781600678212326,
                     0.0081781600678212326, 0.0052165334747187793, 0.0052165334747187793, 0.0022437538722506629,
                     0.0022437538722506629};
      x           = c[n];
    } break;
 
    case 58: {
      T const c[] = {0.0536811198633348488, 0.0536811198633348488, 0.0535263433040582521, 0.0535263433040582521,
                     0.0532172364465790141, 0.0532172364465790141, 0.0527546905263708334, 0.0527546905263708334,
                     0.052140039183669818,  0.052140039183669818,  0.0513750546182857254, 0.0513750546182857254,
                     0.050461942479953125,  0.050461942479953125,  0.049403335508962392,  0.049403335508962392,
                     0.0482022859454177484, 0.0482022859454177484, 0.0468622567290263469, 0.0468622567290263469,
                     0.0453871115148198025, 0.0453871115148198025, 0.0437811035336402510, 0.0437811035336402510,
                     0.0420488633295821259, 0.0420488633295821259, 0.0401953854098677968, 0.0401953854098677968,
                     0.0382260138458584332, 0.0382260138458584332, 0.0361464268670872705, 0.0361464268670872705,
                     0.0339626204934160107, 0.0339626204934160107, 0.0316808912538093273, 0.0316808912538093273,
                     0.0293078180441604907, 0.0293078180441604907, 0.0268502431819818684, 0.0268502431819818684,
                     0.0243152527249639525, 0.0243152527249639525, 0.0217101561401462357, 0.0217101561401462357,
                     0.0190424654618934086, 0.0190424654618934086, 0.0163198742349709650, 0.0163198742349709650,
                     0.0135502371129888121, 0.0135502371129888121, 0.0107415535328787741, 0.0107415535328787741,
                     0.007901973849998674,  0.007901973849998674,  0.0050399816126502430, 0.0050399816126502430,
                     0.0021677232496274499, 0.0021677232496274499};
      x           = c[n];
    } break;
 
    case 59: {
      T const c[] = {0.0527980126219904214, 0.0527244338591279319, 0.0527244338591279319, 0.0525039026478287390,
                     0.0525039026478287390, 0.0521370336483753913, 0.0521370336483753913, 0.0516248493908914821,
                     0.0516248493908914821, 0.0509687774253939168, 0.0509687774253939168, 0.0501706463429969028,
                     0.0501706463429969028, 0.0492326806793619857, 0.0492326806793619857, 0.0481574947146064403,
                     0.0481574947146064403, 0.0469480851869620191, 0.0469480851869620191, 0.045607822940509769,
                     0.045607822940509769,  0.0441404435302973806, 0.0441404435302973806, 0.0425500368110676386,
                     0.0425500368110676386, 0.0408410355386867076, 0.0408410355386867076, 0.039018203016160009,
                     0.039018203016160009,  0.03708661981887092,   0.03708661981887092,   0.0350516696364001087,
                     0.0350516696364001087, 0.0329190242710452777, 0.0329190242710452777, 0.030694627836111683,
                     0.030694627836111683,  0.0283846802005347979, 0.0283846802005347979, 0.0259956197312985001,
                     0.0259956197312985001, 0.023534105393713363,  0.023534105393713363,  0.0210069982884371873,
                     0.0210069982884371873, 0.0184213427536100293, 0.0184213427536100293, 0.0157843473130814661,
                     0.0157843473130814661, 0.0131033663063451910, 0.0131033663063451910, 0.0103858855009958621,
                     0.0103858855009958621, 0.0076395294534875751, 0.0076395294534875751, 0.004872239168265284,
                     0.004872239168265284,  0.0020954922845412234, 0.0020954922845412234};
      x           = c[n];
    } break;
 
    case 60: {
      T const c[] = {0.05190787763122063,   0.05190787763122063,   0.0517679431749101875, 0.0517679431749101875,
                     0.051488451500980933,  0.051488451500980933,  0.0510701560698556274, 0.0510701560698556274,
                     0.0505141845325093745, 0.0505141845325093745, 0.0498220356905501810, 0.0498220356905501810,
                     0.0489955754557568353, 0.0489955754557568353, 0.0480370318199711809, 0.0480370318199711809,
                     0.046948988848912204,  0.046948988848912204,  0.0457343797161144866, 0.0457343797161144866,
                     0.0443964787957871133, 0.0443964787957871133, 0.0429388928359356419, 0.0429388928359356419,
                     0.0413655512355847556, 0.0413655512355847556, 0.0396806954523807994, 0.0396806954523807994,
                     0.0378888675692434440, 0.0378888675692434440, 0.0359948980510845030, 0.0359948980510845030,
                     0.0340038927249464228, 0.0340038927249464228, 0.0319212190192963289, 0.0319212190192963289,
                     0.0297524915007889452, 0.0297524915007889452, 0.0275035567499247916, 0.0275035567499247916,
                     0.025180477621521248,  0.025180477621521248,  0.0227895169439978198, 0.0227895169439978198,
                     0.0203371207294572867, 0.0203371207294572867, 0.0178299010142077202, 0.0178299010142077202,
                     0.0152746185967847993, 0.0152746185967847993, 0.0126781664768159601, 0.0126781664768159601,
                     0.0100475571822879843, 0.0100475571822879843, 0.0073899311633454555, 0.0073899311633454555,
                     0.0047127299269535686, 0.0047127299269535686, 0.0020268119688737584, 0.0020268119688737584};
      x           = c[n];
    } break;
 
    case 61: {
      T const c[] = {0.0510811194407862179, 0.0510144870386972635, 0.0510144870386972635, 0.0508147636688183432,
                     0.0508147636688183432, 0.0504824703867974046, 0.0504824703867974046, 0.0500184741081782534,
                     0.0500184741081782534, 0.0494239853467355899, 0.0494239853467355899, 0.0487005550564115260,
                     0.0487005550564115260, 0.0478500705850956071, 0.0478500705850956071, 0.0468747507508090659,
                     0.0468747507508090659, 0.0457771400531459593, 0.0457771400531459593, 0.0445601020350834882,
                     0.0445601020350834882, 0.0432268118124960979, 0.0432268118124960979, 0.0417807477908884920,
                     0.0417807477908884920, 0.0402256825909982473, 0.0402256825909982473, 0.0385656732070081727,
                     0.0385656732070081727, 0.0368050504231548173, 0.0368050504231548173, 0.0349484075165333510,
                     0.0349484075165333510, 0.0330005882759074106, 0.0330005882759074106, 0.0309666743683973948,
                     0.0309666743683973948, 0.0288519720881834015, 0.0288519720881834015, 0.026661998524150889,
                     0.026661998524150889,  0.0244024671875442029, 0.0244024671875442029, 0.0220792731483190440,
                     0.0220792731483190440, 0.019698477746101181,  0.019698477746101181,  0.0172662929876137435,
                     0.0172662929876137435, 0.0147890658849379145, 0.0147890658849379145, 0.0122732635078121046,
                     0.0122732635078121046, 0.0097254618303561337, 0.0097254618303561337, 0.0071523549917490895,
                     0.0071523549917490895, 0.0045609240060124171, 0.0045609240060124171, 0.0019614533616702826,
                     0.0019614533616702826};
      x           = c[n];
    } break;
 
    case 62: {
      T const c[] = {0.0502480003752562816, 0.0502480003752562816, 0.0501210695690432880, 0.0501210695690432880,
                     0.0498675285949523942, 0.0498675285949523942, 0.0494880179196992925, 0.0494880179196992925,
                     0.0489834962205178371, 0.0489834962205178371, 0.0483552379634776728, 0.0483552379634776728,
                     0.0476048301841012322, 0.0476048301841012322, 0.0467341684784155248, 0.0467341684784155248,
                     0.0457454522145701807, 0.0457454522145701807, 0.0446411789771244142, 0.0446411789771244142,
                     0.0434241382580474195, 0.0434241382580474195, 0.0420974044103850966, 0.0420974044103850966,
                     0.0406643288824174409, 0.0406643288824174409, 0.0391285317519630841, 0.0391285317519630841,
                     0.0374938925822800299, 0.0374938925822800299, 0.0357645406227681412, 0.0357645406227681412,
                     0.0339448443794105450, 0.0339448443794105450, 0.0320394005816246781, 0.0320394005816246781,
                     0.0300530225739898700, 0.0300530225739898700, 0.0279907281633146375, 0.0279907281633146375,
                     0.025857726954024698,  0.025857726954024698,  0.0236594072086827925, 0.0236594072086827925,
                     0.021401322277669968,  0.021401322277669968,  0.0190891766585731987, 0.0190891766585731987,
                     0.0167288117901773162, 0.0167288117901773162, 0.0143261918238065177, 0.0143261918238065177,
                     0.0118873901170105019, 0.0118873901170105019, 0.0094185794284203876, 0.0094185794284203876,
                     0.0069260419018309608, 0.0069260419018309608, 0.0044163334569309048, 0.0044163334569309048,
                     0.0018992056795136904, 0.0018992056795136904};
      x           = c[n];
    } break;
 
    case 63: {
      T const c[] = {0.0494723666239310208, 0.0494118330399181789, 0.0494118330399181789, 0.0492303804237475607,
                     0.0492303804237475607, 0.0489284528205119899, 0.0489284528205119899, 0.0485067890978838478,
                     0.0485067890978838478, 0.0479664211379951314, 0.0479664211379951314, 0.0473086713122689190,
                     0.0473086713122689190, 0.0465351492453836965, 0.0465351492453836965, 0.0456477478762926086,
                     0.0456477478762926086, 0.0446486388259413953, 0.0446486388259413953, 0.0435402670830275907,
                     0.0435402670830275907, 0.042325345020815822,  0.042325345020815822,  0.0410068457596663986,
                     0.0410068457596663986, 0.0395879958915440939, 0.0395879958915440939, 0.0380722675843495567,
                     0.0380722675843495567, 0.0364633700854572896, 0.0364633700854572896, 0.0347652406453558776,
                     0.0347652406453558776, 0.0329820348837793417, 0.0329820348837793417, 0.0311181166222198175,
                     0.0311181166222198175, 0.029178047208280526,  0.029178047208280526,  0.0271665743590979332,
                     0.0271665743590979332, 0.0250886205533449866, 0.0250886205533449866, 0.0229492710048899331,
                     0.0229492710048899331, 0.0207537612580390907, 0.0207537612580390907, 0.018507464460161270,
                     0.018507464460161270,  0.016215878410338338,  0.016215878410338338,  0.0138846126161156108,
                     0.0138846126161156108, 0.011519376076880041,  0.011519376076880041,  0.0091259686763266563,
                     0.0091259686763266563, 0.0067102917659601362, 0.0067102917659601362, 0.0042785083468637618,
                     0.0042785083468637618, 0.001839874595577084,  0.001839874595577084};
      x           = c[n];
    } break;
 
    case 64: {
      T const c[] = {0.0486909570091397203, 0.0486909570091397203, 0.0485754674415034269, 0.0485754674415034269,
                     0.048344762234802957,  0.048344762234802957,  0.0479993885964583077, 0.0479993885964583077,
                     0.0475401657148303086, 0.0475401657148303086, 0.0469681828162100173, 0.0469681828162100173,
                     0.0462847965813144172, 0.0462847965813144172, 0.0454916279274181444, 0.0454916279274181444,
                     0.0445905581637565630, 0.0445905581637565630, 0.0435837245293234533, 0.0435837245293234533,
                     0.0424735151236535890, 0.0424735151236535890, 0.0412625632426235286, 0.0412625632426235286,
                     0.0399537411327203413, 0.0399537411327203413, 0.0385501531786156291, 0.0385501531786156291,
                     0.0370551285402400460, 0.0370551285402400460, 0.0354722132568823838, 0.0354722132568823838,
                     0.0338051618371416093, 0.0338051618371416093, 0.032057928354851553,  0.032057928354851553,
                     0.030234657072402478,  0.030234657072402478,  0.0283396726142594832, 0.0283396726142594832,
                     0.0263774697150546586, 0.0263774697150546586, 0.0243527025687108733, 0.0243527025687108733,
                     0.0222701738083832541, 0.0222701738083832541, 0.0201348231535302093, 0.0201348231535302093,
                     0.0179517157756973430, 0.0179517157756973430, 0.0157260304760247193, 0.0157260304760247193,
                     0.0134630478967186425, 0.0134630478967186425, 0.011168139460131128,  0.011168139460131128,
                     0.0088467598263639477, 0.0088467598263639477, 0.0065044579689783628, 0.0065044579689783628,
                     0.0041470332605624676, 0.0041470332605624676, 0.0017832807216964329, 0.0017832807216964329};
      x           = c[n];
    } break;
  }
  return x;
}
 
} // namespace impl
 
//
// Compute a non-negative integer power of a tensor by binary manipulation.
//
template<typename T, Index N>
KOKKOS_INLINE_FUNCTION
Tensor<T, N>
binary_powering(Tensor<T, N> const & A, Index const exponent)
{
  if (exponent == 0) return eye<T, N>(A.get_dimension());
 
  Index const
  rightmost_bit = 1;
 
  Index const
  number_digits = INDEX_SIZE;
 
  Index const
  leftmost_bit = rightmost_bit << (number_digits - 1);
 
  Index
  t = 0;
 
  for (Index j = 0; j < number_digits; ++j) {
 
    if (((exponent << j) & leftmost_bit) != 0) {
 
      t = number_digits - j - 1;
      break;
 
    }
 
  }
 
  Tensor<T, N>
  P = A;
 
  Index
  i = 0;
 
  Index
  m = exponent;
 
  while ((m & rightmost_bit) == 0) {
    P = P * P;
    ++i;
    m = m >> 1;
  }
 
  Tensor<T, N>
  X = P;
 
  for (Index j = i + 1; j <= t; ++j) {
    P = P * P;
 
    if (((exponent >> j) & rightmost_bit) != 0) {
      X = X * P;
    }
  }
 
  return X;
}
 
//
// Exponential map by squaring and scaling and Padé approximants.
// See algorithm 10.20 in Functions of Matrices, N.J. Higham, SIAM, 2008.
// \param A tensor
// \return \f$ \exp A \f$
//
// Padé truncation-order thresholds of Al-Mohy and Higham, A New Scaling
// and Squaring Algorithm for the Matrix Exponential, SIAM J. Matrix Anal.
// Appl. 31(3), 2009, Table 3.1. The value for order 13 is their revised
// choice that also accounts for rounding in the squaring phase.
template<typename T>
KOKKOS_INLINE_FUNCTION
T
scaling_squaring_theta_ah09(Index const order)
{
  switch (order) {
    default: MT_ERROR_EXIT("Unsupported Pade order for expm."); break;
    case 3: return 1.495585217958292e-02;
    case 5: return 2.539398330063230e-01;
    case 7: return 9.504178996162932e-01;
    case 9: return 2.097847961257068e+00;
    case 13: return 4.25;
  }
  return 0.0;
}
 
// Rounding-error correction term of Al-Mohy and Higham (2009), Eq. (5.1):
// number of extra squarings needed so that the truncation bound remains
// valid when the Padé numerator is evaluated in floating point. Zero for
// all but the most nonnormal arguments.
template<Index N>
Index
expm_ell(Tensor<Real, N> const & A, Index const m)
{
  Index const
  dimension = A.get_dimension();
 
  Real const
  norm_A = norm_1(A);
 
  if (norm_A == 0.0) return 0;
 
  // 1 / |c_{2m+1}| = binomial(2m, m) * (2m+1)!  (few-digit accuracy is
  // all the base-2 logarithm below requires)
  Real
  factorial = 1.0;
 
  for (Index i = 2; i <= 2 * m + 1; ++i) {
    factorial *= static_cast<Real>(i);
  }
 
  Real
  binomial = 1.0;
 
  for (Index i = 1; i <= m; ++i) {
    binomial = binomial * static_cast<Real>(m + i) / static_cast<Real>(i);
  }
 
  Real const
  c_recip = binomial * factorial;
 
  Real const
  u_half = 0.5 * machine_epsilon<Real>();
 
  // conservative bound alpha <= norm(A)^(2m) / c_recip: when it is
  // already below the unit roundoff (in particular whenever norm(A) is
  // within the 2005 norm-based threshold) no correction is needed and
  // the matrix powers below are skipped
  if (std::pow(norm_A, 2.0 * static_cast<Real>(m)) / c_recip <= u_half) {
    return 0;
  }
 
  // norm of |A|^(2m+1), elementwise absolute value, by binary powering
  Tensor<Real, N>
  P(dimension);
 
  for (Index i = 0; i < dimension; ++i) {
    for (Index j = 0; j < dimension; ++j) {
      P(i, j) = std::abs(A(i, j));
    }
  }
 
  Tensor<Real, N> const
  Q = binary_powering(P, 2 * m + 1);
 
  Real const
  alpha = norm_1(Q) / (norm_A * c_recip);
 
  if (alpha == 0.0) return 0;
 
  int const
  extra = static_cast<int>(
      std::ceil(std::log2(alpha / u_half) / (2.0 * static_cast<Real>(m))));
 
  return extra > 0 ? static_cast<Index>(extra) : 0;
}
 
template <typename T, Index N> Tensor<T, N> exp_pade(Tensor<T, N> const &A) {
  // Scaling and squaring with diagonal Padé approximants and the
  // parameter selection of Al-Mohy and Higham (2009). Their theta_m
  // thresholds satisfy theta_m^(2m) / c_recip(m) <= u, so whenever
  // norm(A) <= theta_m both the eta_m test (eta_m <= norm) and the
  // vanishing of the expm_ell correction are guaranteed: for such
  // arguments selecting the order directly from norm(A) is exact
  // short-circuit evaluation of the full algorithm at a fraction of the
  // cost. The norm(A^p)^(1/p) analysis, which prevents overscaling of
  // nonnormal arguments, runs only when norm(A) > theta_13.
  Index const
  dimension = A.get_dimension();
 
  Real const
  norm = Sacado::ScalarValue<T>::eval(norm_1(A));
 
  Index
  order = 13;
 
  Index const
  orders[] = {3, 5, 7, 9, 13};
 
  for (Index i = 0; i < 5; ++i) {
    if (norm <= scaling_squaring_theta_ah09<Real>(orders[i])) {
      order = orders[i];
      break;
    }
  }
 
  Real
  eta_5 = 0.0;
 
  bool const
  scaling_needed = norm > scaling_squaring_theta_ah09<Real>(13);
 
  bool
  use_ell = false;
 
  Tensor<Real, N>
  Ar(dimension);
 
  if (scaling_needed) {
    // Double-precision value copy for parameter selection only; carries
    // no derivative information.
    for (Index i = 0; i < dimension; ++i) {
      for (Index j = 0; j < dimension; ++j) {
        Ar(i, j) = Sacado::ScalarValue<T>::eval(A(i, j));
      }
    }
    Tensor<Real, N> const A2r = Ar * Ar;
    Tensor<Real, N> const A4r = A2r * A2r;
    Tensor<Real, N> const A6r = A2r * A4r;
    Tensor<Real, N> const A8r = A4r * A4r;
    Real const d6 = std::pow(norm_1(A6r), 1.0 / 6.0);
    Real const d8 = std::pow(norm_1(A8r), 1.0 / 8.0);
    Real const d10 = std::pow(norm_1(A4r * A6r), 1.0 / 10.0);
    Real const eta_3 = std::max(d6, d8);
    Real const eta_4 = std::max(d8, d10);
    eta_5 = std::min(eta_3, eta_4);
    use_ell = true;
  }
 
  Tensor<T, N>
  B;
 
  if (order < 13) {
 
    Tensor<T, N>
    U;
 
    Tensor<T, N>
    V;
 
    std::tie(U, V) = impl::pade_polynomial_terms(A, order);
 
    // The Pade denominator V - U = q_m(A) is diagonally dominant and
    // well conditioned for norm(A) within the theta thresholds, so the
    // inexpensive adjugate inverse is safe here.
    B = inverse(V - U) * (U + V);
 
  } else {
 
      Real const
      theta_highest = scaling_squaring_theta_ah09<Real>(13);
 
      int
      signed_power = 0;
 
      if (eta_5 > 0.0) {
        signed_power = static_cast<int>(
            std::ceil(std::log2(eta_5 / theta_highest)));
        if (signed_power < 0) signed_power = 0;
      }
 
      // rounding correction on the scaled argument (needed only when the
      // alpha-based scaling was used; the norm-based bounds make it zero)
      if (use_ell) {
        Real
        scale_r = 1.0;
        for (int j = 0; j < signed_power; ++j) {
          scale_r /= 2.0;
        }
        signed_power += static_cast<int>(expm_ell(scale_r * Ar, 13));
      }
 
      Index const
      power_two = signed_power > 0 ? static_cast<Index>(signed_power) : 0;
 
      Real
      scale = 1.0;
 
      for (Index j = 0; j < power_two; ++j) {
        scale /= 2.0;
      }
 
      Tensor<T, N> const
      I = identity<T, N>(dimension);
 
      Tensor<T, N> const
      A1 = scale * A;
 
      Tensor<T, N> const
      A2 = A1 * A1;
 
      Tensor<T, N> const
      A4 = A2 * A2;
 
      Tensor<T, N> const
      A6 = A2 * A4;
 
      Real const b0  = impl::polynomial_coefficient<Real>(order, 0);
      Real const b1  = impl::polynomial_coefficient<Real>(order, 1);
      Real const b2  = impl::polynomial_coefficient<Real>(order, 2);
      Real const b3  = impl::polynomial_coefficient<Real>(order, 3);
      Real const b4  = impl::polynomial_coefficient<Real>(order, 4);
      Real const b5  = impl::polynomial_coefficient<Real>(order, 5);
      Real const b6  = impl::polynomial_coefficient<Real>(order, 6);
      Real const b7  = impl::polynomial_coefficient<Real>(order, 7);
      Real const b8  = impl::polynomial_coefficient<Real>(order, 8);
      Real const b9  = impl::polynomial_coefficient<Real>(order, 9);
      Real const b10 = impl::polynomial_coefficient<Real>(order, 10);
      Real const b11 = impl::polynomial_coefficient<Real>(order, 11);
      Real const b12 = impl::polynomial_coefficient<Real>(order, 12);
      Real const b13 = impl::polynomial_coefficient<Real>(order, 13);
 
      Tensor<T, N> const
      U = A1 * (
          (A6 * (b13 * A6 + b11 * A4 + b9 * A2) +
              b7 * A6 + b5 * A4 + b3 * A2 + b1 * I));
 
      Tensor<T, N> const
      V = A6 * (b12 * A6 + b10 * A4 + b8 * A2) +
      b6 * A6 + b4 * A4 + b2 * A2 + b0 * I;
 
      Tensor<T, N> const
      R = inverse(V - U) * (U + V);
 
      Index const
      exponent = (1U << power_two);
 
      B = binary_powering(R, exponent);
 
  }
 
  return B;
}
 
//
// Logarithmic map by Taylor series.
//
template<typename T, Index N>
KOKKOS_INLINE_FUNCTION
Tensor<T, N>
log_taylor(Tensor<T, N> const & A)
{
  Index const
  max_iter = 128;
 
  T const
  tol = machine_epsilon<T>();
 
  T const
  norm_tensor = norm_1(A);
 
  Index const
  dimension = A.get_dimension();
 
  Tensor<T, N> const
  A_minus_I = A - identity<T, N>(dimension);
 
  Tensor<T, N>
  term = A_minus_I;
 
  T
  norm_term = norm_1(term);
 
  T
  relative_error = norm_term / norm_tensor;
 
  Tensor<T, N>
  B = term;
 
  Index
  k = 1;
 
  while (relative_error > tol && k <= max_iter) {
    term = static_cast<T>(- (k / (k + 1.0))) * term * A_minus_I;
    B = B + term;
    norm_term = norm_1(term);
    relative_error = norm_term / norm_tensor;
    ++k;
  }
 
  return B;
}
 
//
// Logarithmic map.
//
template<typename T, Index N>
KOKKOS_INLINE_FUNCTION
Tensor<T, N>
log(Tensor<T, N> const & A)
{
  return log_pade(A);
}
 
//
// Logarithmic map by Gregory series.
//
template<typename T, Index N>
KOKKOS_INLINE_FUNCTION
Tensor<T, N>
log_gregory(Tensor<T, N> const & A)
{
  Index const
  max_iter = 8192;
 
  T const
  tol = machine_epsilon<T>();
 
  T const
  norm_tensor = norm_1(A);
 
  Index const
  dimension = A.get_dimension();
 
  Tensor<T, N> const
  I_minus_A = identity<T, N>(dimension) - A;
 
  Tensor<T, N> const
  I_plus_A = identity<T, N>(dimension) + A;
 
  Tensor<T, N>
  term = I_minus_A * inverse(I_plus_A);
 
  T
  norm_term = norm_1(term);
 
  T
  relative_error = norm_term / norm_tensor;
 
  Tensor<T, N> const
  C = term * term;
 
  Tensor<T, N>
  B = term;
 
  Index
  k = 1;
 
  while (relative_error > tol && k <= max_iter + 1) {
    term = static_cast<T>((2 * k - 1.0) / (2 * k + 1.0)) * term * C;
    B = B + term;
    norm_term = norm_1(term);
    relative_error = norm_term / norm_tensor;
    ++k;
  }
 
  B = - 2.0 * B;
 
  return B;
}
 
// Matrix square root by product form of Denman-Beavers iteration.
template<typename T, Index N>
KOKKOS_INLINE_FUNCTION
Tensor<T, N>
sqrt_dbp(Tensor<T, N> const & A, int& k)
{
  auto const dimension = A.get_dimension();
  auto const eps = machine_epsilon<T>();
  auto const tol = 0.5 * std::sqrt(dimension) * eps;
  auto const I = identity<T, N>(dimension);
  auto const max_iter = 32;
  auto X = A;
  auto M = A;
  auto scale = true;
  k = 0;
  while (k++ < max_iter) {
    if (scale == true) {
      auto const d  = abs(det(M));
      auto const d2 = std::sqrt(d);
      auto const d6 = cbrt(d2);
      auto const g  = 1.0 / d6;
      // For ill-conditioned M the expanded determinant can cancel to zero
      // (or overflow), which would make g non-finite and flood the iteration
      // with Inf/NaN. The scaling only accelerates convergence, so it is
      // safe to skip it in that case.
      if (Sacado::ScalarValue<T>::eval(d) > 0.0 &&
          std::isfinite(Sacado::ScalarValue<T>::eval(g))) {
        X *= g;
        M *= g * g;
      }
    }
    auto const Y = X;
    auto const L = inverse_full_pivot(M);
    X = 0.5 * X * (I + L);
    M = 0.5 * (I + 0.5 * (M + L));
    auto const error = norm(M - I);
    auto const diff  = norm(X - Y) / norm(X);
    scale            = diff >= 0.01;
    if (error <= tol) break;
  }
  return X;
}
 
// Tensor square root
template<typename T, Index N>
KOKKOS_INLINE_FUNCTION
Tensor<T, N>
sqrt(Tensor<T, N> const & A)
{
  int k = 0;
  return sqrt_dbp(A, k);
}
 
// Logarithmic map by Padé approximant and partial fractions
template<typename T, Index N>
KOKKOS_INLINE_FUNCTION
Tensor<T, N>
log_pade_pf(Tensor<T, N> const & A, Index const n)
{
  auto const dimension = A.get_dimension();
  auto const I = identity<T, N>(dimension);
  auto       X = 0.0 * A;
  for (Index i = 0; i < n; ++i) {
    auto const x = 0.5 * (1.0 + impl::gauss_legendre_abscissae<T>(n, i));
    auto const w = 0.5 * impl::gauss_legendre_weights<T>(n, i);
    auto const B = I + x * A;
    X += w * A * inverse_full_pivot(B);
  }
  return X;
}
 
// Logarithmic map by inverse scaling and squaring and Padé approximants
template<typename T, Index N>
KOKKOS_INLINE_FUNCTION
Tensor<T, N>
log_iss(Tensor<T, N> const & A)
{
  auto const dimension = A.get_dimension();
  auto const I   = identity<T, N>(dimension);
  auto const c15 = impl::pade_coefficients<T>(15);
  auto       X   = A;
  auto       i   = 5;
  auto       j   = 0;
  auto       k   = 0;
  auto       m   = 0;
  while (true) {
    auto const diff = norm_1(X - I);
    if (diff <= c15) {
      auto p = 2;
      while (impl::pade_coefficients<T>(p) <= diff && p < 16) {
        ++p;
      }
      auto q = 2;
      while (impl::pade_coefficients<T>(q) <= diff / 2.0 && q < 16) {
        ++q;
      }
      if ((2 * (p - q) / 3) < i || ++j == 2) {
        m = p + 1;
        break;
      }
    }
    X = sqrt_dbp(X, i);
    ++k;
  }
  X = (1U << k) * log_pade_pf(X - I, m);
  return X;
}
 
namespace impl {
 
// Givens pair (c, s) with G = [c, -s; s, c] such that G [a; b] = [r; 0].
template<typename T>
KOKKOS_INLINE_FUNCTION
std::pair<T, T>
givens_zero(T const & a, T const & b)
{
  T const r = std::sqrt(a * a + b * b);
  if (Sacado::ScalarValue<T>::eval(r) == 0.0) {
    return std::make_pair(T(1.0), T(0.0));
  }
  return std::make_pair(a / r, -b / r);
}
 
// Standardize the 2x2 diagonal block of S at rows/columns (p, p+1) by an
// orthogonal similarity applied to all of S and accumulated into Q:
// afterwards the block is upper triangular (real eigenvalues, S(p+1,p)
// set exactly to zero) or has equal diagonal entries and off-diagonal
// entries of opposite sign (complex conjugate pair), as in LAPACK dlanv2.
template<typename T, Index N>
void
schur_standardize_2x2(Tensor<T, N> & S, Tensor<T, N> & Q, Index const p)
{
  Index const q = p + 1;
 
  T const a = S(p, p);
  T const b = S(p, q);
  T const c = S(q, p);
  T const d = S(q, q);
 
  if (Sacado::ScalarValue<T>::eval(std::abs(c)) == 0.0) {
    S(q, p) = 0.0;
    return;
  }
 
  T const half_diff = 0.5 * (a - d);
  T const disc = half_diff * half_diff + b * c;
 
  T cs, sn;
 
  if (Sacado::ScalarValue<T>::eval(disc) >= 0.0) {
    // Real eigenvalues: rotate the eigenvector of the larger-|z| root
    // onto the first coordinate axis.
    T const root = std::sqrt(disc);
    T const z = Sacado::ScalarValue<T>::eval(half_diff) >= 0.0 ?
        T(half_diff + root) : T(half_diff - root);
    std::tie(cs, sn) = givens_zero(z, c);
    givens_left(cs, sn, p, q, S);
    givens_right(cs, sn, p, q, S);
    givens_right(cs, sn, p, q, Q);
    S(q, p) = 0.0;
  } else {
    // Complex pair: rotate so the diagonal entries are equal.
    T const theta = 0.5 * std::atan2(a - d, b + c);
    cs = std::cos(theta);
    sn = std::sin(theta);
    givens_left(cs, sn, p, q, S);
    givens_right(cs, sn, p, q, S);
    givens_right(cs, sn, p, q, Q);
    T const m = 0.5 * (S(p, p) + S(q, q));
    S(p, p) = m;
    S(q, q) = m;
  }
}
 
// Real Schur decomposition A = Q S Q^T for dimension <= 3: S is upper
// quasi-triangular with 1x1 blocks and standardized 2x2 blocks (complex
// conjugate eigenvalue pairs), Q is orthogonal. Hessenberg reduction by a
// Givens rotation followed by Francis double-shift QR iteration.
template<typename T, Index N>
std::pair<Tensor<T, N>, Tensor<T, N>>
schur_real_small(Tensor<T, N> const & A)
{
  Index const
  dimension = A.get_dimension();
 
  Tensor<T, N>
  S = A;
 
  Tensor<T, N>
  Q = identity<T, N>(dimension);
 
  if (dimension == 1) {
    return std::make_pair(Q, S);
  }
 
  if (dimension == 2) {
    schur_standardize_2x2(S, Q, 0);
    return std::make_pair(Q, S);
  }
 
  // dimension == 3
 
  // Hessenberg: zero S(2,0) with a rotation in the (1,2) plane.
  {
    auto const cs_sn = givens_zero(S(1, 0), S(2, 0));
    givens_left(cs_sn.first, cs_sn.second, 1, 2, S);
    givens_right(cs_sn.first, cs_sn.second, 1, 2, S);
    givens_right(cs_sn.first, cs_sn.second, 1, 2, Q);
    S(2, 0) = 0.0;
  }
 
  Real const
  eps = machine_epsilon<Real>();
 
  Index const
  max_iter = 64;
 
  for (Index iter = 0; iter < max_iter; ++iter) {
 
    Real const h00 = std::abs(Sacado::ScalarValue<T>::eval(S(0, 0)));
    Real const h11 = std::abs(Sacado::ScalarValue<T>::eval(S(1, 1)));
    Real const h22 = std::abs(Sacado::ScalarValue<T>::eval(S(2, 2)));
    Real const h10 = std::abs(Sacado::ScalarValue<T>::eval(S(1, 0)));
    Real const h21 = std::abs(Sacado::ScalarValue<T>::eval(S(2, 1)));
 
    if (h10 <= eps * (h00 + h11)) {
      // leading 1x1 deflated; standardize the trailing 2x2
      S(1, 0) = 0.0;
      schur_standardize_2x2(S, Q, 1);
      return std::make_pair(Q, S);
    }
 
    if (h21 <= eps * (h11 + h22)) {
      // trailing 1x1 deflated; standardize the leading 2x2
      S(2, 1) = 0.0;
      schur_standardize_2x2(S, Q, 0);
      return std::make_pair(Q, S);
    }
 
    // Francis double-shift step: first column of (S - l1 I)(S - l2 I)
    // with l1, l2 the eigenvalues of the trailing 2x2 block, written in
    // the classical difference form so that clustered eigenvalues (for
    // example A near the identity, where the column is O(eps^2) but the
    // naive expansion sums O(1) terms) do not lose the shift direction to
    // cancellation. Every 12th iteration perturb the shift to escape
    // stagnation.
    T x, y, z;
 
    if (iter > 0 && iter % 12 == 0) {
      T const shift = std::abs(S(2, 1)) + std::abs(S(1, 0));
      T const tr = 2.0 * shift;
      T const dt = shift * shift;
      x = S(0, 0) * S(0, 0) + S(0, 1) * S(1, 0) - tr * S(0, 0) + dt;
      y = S(1, 0) * (S(0, 0) + S(1, 1) - tr);
      z = S(1, 0) * S(2, 1);
    } else {
      x = (S(0, 0) - S(1, 1)) * (S(0, 0) - S(2, 2)) -
          S(1, 2) * S(2, 1) + S(0, 1) * S(1, 0);
      y = S(1, 0) * (S(0, 0) - S(2, 2));
      z = S(1, 0) * S(2, 1);
    }
 
    // Two rotations mapping (x, y, z) onto e1
    auto const g1 = givens_zero(y, z);
    T const y1 = g1.first * y - g1.second * z;
    auto const g2 = givens_zero(x, y1);
 
    givens_left(g1.first, g1.second, 1, 2, S);
    givens_right(g1.first, g1.second, 1, 2, S);
    givens_right(g1.first, g1.second, 1, 2, Q);
 
    givens_left(g2.first, g2.second, 0, 1, S);
    givens_right(g2.first, g2.second, 0, 1, S);
    givens_right(g2.first, g2.second, 0, 1, Q);
 
    // Restore Hessenberg form: zero the bulge S(2,0)
    auto const g3 = givens_zero(S(1, 0), S(2, 0));
    givens_left(g3.first, g3.second, 1, 2, S);
    givens_right(g3.first, g3.second, 1, 2, S);
    givens_right(g3.first, g3.second, 1, 2, Q);
    S(2, 0) = 0.0;
  }
 
  MT_ERROR_EXIT("Schur QR iteration failed to converge.");
  return std::make_pair(Q, S);
}
 
// Principal square root of an upper quasi-triangular S (dimension <= 3)
// with standardized 2x2 blocks and spectrum off the closed negative real
// axis. 1x1 blocks: scalar square root. 2x2 blocks [m, b; c, m], bc < 0:
// closed form sqrt(B) = (B + delta I) / tau with delta = sqrt(det B) and
// tau = sqrt(tr B + 2 delta). Off-diagonal blocks from the Sylvester
// equations of U^2 = S.
template<typename T, Index N>
Tensor<T, N>
sqrt_quasi_triu(Tensor<T, N> const & S)
{
  Index const
  dimension = S.get_dimension();
 
  Tensor<T, N>
  U(dimension, Filler::ZEROS);
 
  // block structure from exact subdiagonal zeros
  Index
  blocks[3];
 
  Index
  starts[3];
 
  Index
  num_blocks = 0;
 
  for (Index i = 0; i < dimension;) {
    bool const two = (i + 1 < dimension) &&
        (Sacado::ScalarValue<T>::eval(S(i + 1, i)) != 0.0);
    starts[num_blocks] = i;
    blocks[num_blocks] = two ? 2 : 1;
    ++num_blocks;
    i += two ? 2 : 1;
  }
 
  // diagonal blocks
  for (Index b = 0; b < num_blocks; ++b) {
    Index const i = starts[b];
    if (blocks[b] == 1) {
      U(i, i) = std::sqrt(S(i, i));
    } else {
      T const m = S(i, i);
      T const delta =
          std::sqrt(m * m - S(i, i + 1) * S(i + 1, i));
      T const tau = std::sqrt(m + m + delta + delta);
      U(i, i) = (m + delta) / tau;
      U(i + 1, i + 1) = U(i, i);
      U(i, i + 1) = S(i, i + 1) / tau;
      U(i + 1, i) = S(i + 1, i) / tau;
    }
  }
 
  // off-diagonal blocks, nearest first
  for (Index gap = 1; gap < num_blocks; ++gap) {
    for (Index b = 0; b + gap < num_blocks; ++b) {
      Index const bi = b;
      Index const bj = b + gap;
      Index const i = starts[bi];
      Index const j = starts[bj];
      // right-hand side: S_ij - sum_k U_ik U_kj over blocks strictly
      // between bi and bj
      T rhs[2][2] = {{T(0.0), T(0.0)}, {T(0.0), T(0.0)}};
      for (Index r = 0; r < blocks[bi]; ++r) {
        for (Index c = 0; c < blocks[bj]; ++c) {
          rhs[r][c] = S(i + r, j + c);
        }
      }
      for (Index bk = bi + 1; bk < bj; ++bk) {
        Index const k = starts[bk];
        for (Index r = 0; r < blocks[bi]; ++r) {
          for (Index c = 0; c < blocks[bj]; ++c) {
            for (Index kk = 0; kk < blocks[bk]; ++kk) {
              rhs[r][c] -= U(i + r, k + kk) * U(k + kk, j + c);
            }
          }
        }
      }
      // solve U_ii X + X U_jj = rhs
      if (blocks[bi] == 1 && blocks[bj] == 1) {
        U(i, j) = rhs[0][0] / (U(i, i) + U(j, j));
      } else if (blocks[bi] == 2 && blocks[bj] == 1) {
        // (U_ii + u_jj I) x = rhs (2x2 solve)
        T const a00 = U(i, i) + U(j, j);
        T const a01 = U(i, i + 1);
        T const a10 = U(i + 1, i);
        T const a11 = U(i + 1, i + 1) + U(j, j);
        T const det = a00 * a11 - a01 * a10;
        U(i, j) = (rhs[0][0] * a11 - a01 * rhs[1][0]) / det;
        U(i + 1, j) = (a00 * rhs[1][0] - a10 * rhs[0][0]) / det;
      } else {
        // x (U_jj + u_ii I) = rhs (1x2 row solve)
        T const a00 = U(j, j) + U(i, i);
        T const a01 = U(j, j + 1);
        T const a10 = U(j + 1, j);
        T const a11 = U(j + 1, j + 1) + U(i, i);
        T const det = a00 * a11 - a01 * a10;
        U(i, j) = (rhs[0][0] * a11 - rhs[0][1] * a10) / det;
        U(i, j + 1) = (rhs[0][1] * a00 - rhs[0][0] * a01) / det;
      }
    }
  }
 
  return U;
}
 
} // namespace impl
 
//
// Logarithmic map by real Schur decomposition and inverse scaling and
// squaring on the quasi-triangular factor, in the spirit of Al-Mohy and
// Higham, Improved Inverse Scaling and Squaring Algorithms for the
// Matrix Logarithm, SIAM J. Sci. Comput. 34(4), 2012. Dimension <= 3.
//
template<typename T, Index N>
Tensor<T, N>
log_schur(Tensor<T, N> const & A)
{
  Index const
  dimension = A.get_dimension();
 
  Tensor<T, N>
  Q(dimension);
 
  Tensor<T, N>
  R(dimension);
 
  std::tie(Q, R) = impl::schur_real_small(A);
 
  // The principal logarithm requires the spectrum off the closed
  // negative real axis; real eigenvalues appear as 1x1 diagonal blocks.
  for (Index i = 0; i < dimension; ++i) {
    bool const in_block =
        (i + 1 < dimension &&
         Sacado::ScalarValue<T>::eval(R(i + 1, i)) != 0.0) ||
        (i > 0 && Sacado::ScalarValue<T>::eval(R(i, i - 1)) != 0.0);
    if (!in_block && Sacado::ScalarValue<T>::eval(R(i, i)) <= 0.0) {
      MT_ERROR_EXIT("Non-positive real eigenvalue: principal log undefined.");
    }
  }
 
  // Inverse scaling and squaring on R with exact quasi-triangular square
  // roots; degree selection as in log_iss.
  auto const
  I = identity<T, N>(dimension);
 
  auto const
  c15 = impl::pade_coefficients<Real>(15);
 
  Tensor<T, N>
  X = R;
 
  auto j = 0;
  auto k = 0;
  auto m = 0;
 
  while (true) {
    Real const
    diff = Sacado::ScalarValue<T>::eval(norm_1(X - I));
    if (diff <= c15) {
      auto p = 2;
      while (impl::pade_coefficients<Real>(p) <= diff && p < 16) {
        ++p;
      }
      auto q = 2;
      while (impl::pade_coefficients<Real>(q) <= diff / 2.0 && q < 16) {
        ++q;
      }
      // an exact triangular square root is cheap: take another one only
      // if it reduces the Pade degree meaningfully
      if (p - q < 3 || ++j == 2) {
        m = p + 1;
        break;
      }
    }
    X = impl::sqrt_quasi_triu(X);
    ++k;
  }
 
  Tensor<T, N>
  L = (1U << k) * log_pade_pf(X - I, m);
 
  // Recompute the diagonal blocks of log R exactly from the eigenvalues
  // (Al-Mohy-Higham refinement): 1x1 blocks are scalar logs; for a
  // standardized 2x2 block [m, b; c, m] with eigenvalues m +/- i mu,
  // mu = sqrt(-b c), the log is l I + (theta / mu) K with
  // l = log(sqrt(m^2 + mu^2)), theta = atan2(mu, m), K = [0, b; c, 0].
  for (Index i = 0; i < dimension;) {
    bool const two = (i + 1 < dimension) &&
        (Sacado::ScalarValue<T>::eval(R(i + 1, i)) != 0.0);
    if (!two) {
      L(i, i) = std::log(R(i, i));
      i += 1;
    } else {
      T const mm = R(i, i);
      T const mu = std::sqrt(-R(i, i + 1) * R(i + 1, i));
      T const ell = 0.5 * std::log(mm * mm + mu * mu);
      T const theta = std::atan2(mu, mm);
      T const factor = theta / mu;
      L(i, i) = ell;
      L(i + 1, i + 1) = ell;
      L(i, i + 1) = R(i, i + 1) * factor;
      L(i + 1, i) = R(i + 1, i) * factor;
      i += 2;
    }
  }
 
  return dot(Q, dot_t(L, Q));
}
 
//
// Logarithmic map by squaring and scaling and Padé approximants.
// See algorithm 11.10 in Functions of Matrices, N.J. Higham, SIAM, 2008.
// \param A tensor
// \return \f$ \log A \f$
//
template<typename T, Index N>
KOKKOS_INLINE_FUNCTION
Tensor<T, N>
log_pade(Tensor<T, N> const & A)
{
  Index const
  dimension = A.get_dimension();
 
  if (dimension <= 3) {
 
    // Near the identity no square roots are needed: evaluate the Padé
    // approximant directly on A - I, which is computed exactly. This
    // path preserves the relative accuracy of small logarithms, which
    // any similarity-based method loses to the ~eps*norm(A) rounding of
    // the transformation itself.
    Tensor<T, N> const
    W = A - identity<T, N>(dimension);
 
    Real const
    diff = Sacado::ScalarValue<T>::eval(norm_1(W));
 
    if (diff <= impl::pade_coefficients<Real>(15)) {
      auto p = 2;
      while (impl::pade_coefficients<Real>(p) <= diff && p < 16) {
        ++p;
      }
      return log_pade_pf(W, p + 1);
    }
 
    // Away from the identity: Schur-based inverse scaling and squaring.
    return log_schur(A);
  }
 
  return log_iss(A);
}
 
//
// Logarithmic map for symmetric tensor.
//
template<typename T, Index N>
KOKKOS_INLINE_FUNCTION
Tensor<T, N>
log_sym(Tensor<T, N> const & A)
{
  return log_eig_sym(A);
}
 
//
// Logarithmic map for symmetric tensor using eigenvalue decomposition.
//
template<typename T, Index N>
KOKKOS_INLINE_FUNCTION
Tensor<T, N>
log_eig_sym(Tensor<T, N> const & A)
{
  Index const
  dimension = A.get_dimension();
 
  Tensor<T, N>
  V(dimension);
 
  Tensor<T, N>
  D(dimension);
 
  std::tie(V, D) = eig_sym(A);
 
  for (Index i = 0; i < dimension; ++i) {
    if (D(i, i) < T(0)) {
      MT_ERROR_EXIT("Non-SPD input: negative eigenvalue.");
    }
    D(i, i) = std::log(D(i, i));
  }
 
  Tensor<T, N> const
  B = dot_t(dot(V, D), V);
 
  return B;
}
 
//
// Exponential map for symmetric tensor.
//
template<typename T, Index N>
KOKKOS_INLINE_FUNCTION
Tensor<T, N>
exp_sym(Tensor<T, N> const & A)
{
  return exp_eig_sym(A);
}
 
//
// Exponential map for symmetric tensor using eigenvalue decomposition.
//
template<typename T, Index N>
KOKKOS_INLINE_FUNCTION
Tensor<T, N>
exp_eig_sym(Tensor<T, N> const & A)
{
  Index const
  dimension = A.get_dimension();
 
  Tensor<T, N>
  V(dimension);
 
  Tensor<T, N>
  D(dimension);
 
  std::tie(V, D) = eig_sym(A);
 
  for (Index i = 0; i < dimension; ++i) {
    D(i, i) = std::exp(D(i, i));
  }
 
  Tensor<T, N> const
  B = dot_t(dot(V, D), V);
 
  return B;
}
 
//
// R^N logarithmic map using BCH expansion (4 terms)
// \param x tensor
// \param y tensor
// \return Baker-Campbell-Hausdorff series up to 4 terms
//
template<typename T, Index N>
KOKKOS_INLINE_FUNCTION
Tensor<T, N>
bch(Tensor<T, N> const & x, Tensor<T, N> const & y)
{
  return
      // first order term
      x + y
      +
      // second order term
      0.5*(x*y - y*x)
      +
      // third order term
      1.0/12.0 *
      (x*x*y - 2.0*x*y*x + x*y*y + y*x*x - 2.0*y*x*y + y*y*x)
      +
      // fourth order term
      1.0/24.0 *
      (x*x*y*y - 2.0*x*y*x*y + 2.0*y*x*y*x - y*y*x*x);
}
 
} // namespace minitensor
 
#endif // MiniTensor_MatrixFunctions_h