1 // This code has been mechanically translated from the original FORTRAN 2 // code at http://netlib.org/quadpack. 3 4 /** Authors: Lars Tandle Kyllingstad 5 Copyright: Copyright (c) 2009, Lars T. Kyllingstad. All rights reserved. 6 License: Boost License 1.0 7 */ 8 module scid.ports.quadpack.qcheb; 9 10 11 import scid.core.fortran; 12 13 14 15 16 /// 17 void qcheb(Real)(const Real* x_, Real* fval_, Real* cheb12_, Real* cheb24_) 18 { 19 //***begin prologue dqcheb 20 //***refer to dqc25c,dqc25f,dqc25s 21 //***routines called (none) 22 //***revision date 830518 (yymmdd) 23 //***keywords chebyshev series expansion, fast fourier transform 24 //***author piessens,robert,appl. math. & progr. div. - k.u.leuven 25 // de doncker,elise,appl. math. & progr. div. - k.u.leuven 26 //***purpose this routine computes the chebyshev series expansion 27 // of degrees 12 and 24 of a function using a 28 // fast fourier transform method 29 // f(x) = sum(k=1,..,13) (cheb12(k)*t(k-1,x)), 30 // f(x) = sum(k=1,..,25) (cheb24(k)*t(k-1,x)), 31 // where t(k,x) is the chebyshev polynomial of degree k. 32 //***description 33 // 34 // chebyshev series expansion 35 // standard fortran subroutine 36 // double precision version 37 // 38 // parameters 39 // on entry 40 // x - double precision 41 // vector of dimension 11 containing the 42 // values cos(k*pi/24), k = 1, ..., 11 43 // 44 // fval - double precision 45 // vector of dimension 25 containing the 46 // function values at the points 47 // (b+a+(b-a)*cos(k*pi/24))/2, k = 0, ...,24, 48 // where (a,b) is the approximation interval. 49 // fval(1) and fval(25) are divided by two 50 // (these values are destroyed at output). 51 // 52 // on return 53 // cheb12 - double precision 54 // vector of dimension 13 containing the 55 // chebyshev coefficients for degree 12 56 // 57 // cheb24 - double precision 58 // vector of dimension 25 containing the 59 // chebyshev coefficients for degree 24 60 // 61 //***end prologue dqcheb 62 // 63 Real alam,alam1,alam2,part1,part2, part3; 64 Real[12] v_; 65 int i=1,j=1; 66 // 67 auto cheb12 = dimension(cheb12_, 13); 68 auto cheb24 = dimension(cheb24_, 25); 69 auto fval = dimension(fval_, 25); 70 auto v = dimension(v_.ptr, 12); 71 auto x = dimension(x_, 11); 72 // 73 //***first executable statement dqcheb 74 for (i=1; i<=12; i++) { // end: 10 75 j = 26-i; 76 v[i] = fval[i]-fval[j]; 77 fval[i] = fval[i]+fval[j]; 78 l10: ;} 79 alam1 = v[1]-v[9]; 80 alam2 = x[6]*(v[3]-v[7]-v[11]); 81 cheb12[4] = alam1+alam2; 82 cheb12[10] = alam1-alam2; 83 alam1 = v[2]-v[8]-v[10]; 84 alam2 = v[4]-v[6]-v[12]; 85 alam = x[3]*alam1+x[9]*alam2; 86 cheb24[4] = cheb12[4]+alam; 87 cheb24[22] = cheb12[4]-alam; 88 alam = x[9]*alam1-x[3]*alam2; 89 cheb24[10] = cheb12[10]+alam; 90 cheb24[16] = cheb12[10]-alam; 91 part1 = x[4]*v[5]; 92 part2 = x[8]*v[9]; 93 part3 = x[6]*v[7]; 94 alam1 = v[1]+part1+part2; 95 alam2 = x[2]*v[3]+part3+x[10]*v[11]; 96 cheb12[2] = alam1+alam2; 97 cheb12[12] = alam1-alam2; 98 alam = x[1]*v[2]+x[3]*v[4]+x[5]*v[6]+x[7]*v[8] 99 +x[9]*v[10]+x[11]*v[12]; 100 cheb24[2] = cheb12[2]+alam; 101 cheb24[24] = cheb12[2]-alam; 102 alam = x[11]*v[2]-x[9]*v[4]+x[7]*v[6]-x[5]*v[8] 103 +x[3]*v[10]-x[1]*v[12]; 104 cheb24[12] = cheb12[12]+alam; 105 cheb24[14] = cheb12[12]-alam; 106 alam1 = v[1]-part1+part2; 107 alam2 = x[10]*v[3]-part3+x[2]*v[11]; 108 cheb12[6] = alam1+alam2; 109 cheb12[8] = alam1-alam2; 110 alam = x[5]*v[2]-x[9]*v[4]-x[1]*v[6] 111 -x[11]*v[8]+x[3]*v[10]+x[7]*v[12]; 112 cheb24[6] = cheb12[6]+alam; 113 cheb24[20] = cheb12[6]-alam; 114 alam = x[7]*v[2]-x[3]*v[4]-x[11]*v[6]+x[1]*v[8] 115 -x[9]*v[10]-x[5]*v[12]; 116 cheb24[8] = cheb12[8]+alam; 117 cheb24[18] = cheb12[8]-alam; 118 for (i=1; i<=6; i++) { // end: 20 119 j = 14-i; 120 v[i] = fval[i]-fval[j]; 121 fval[i] = fval[i]+fval[j]; 122 l20: ;} 123 alam1 = v[1]+x[8]*v[5]; 124 alam2 = x[4]*v[3]; 125 cheb12[3] = alam1+alam2; 126 cheb12[11] = alam1-alam2; 127 cheb12[7] = v[1]-v[5]; 128 alam = x[2]*v[2]+x[6]*v[4]+x[10]*v[6]; 129 cheb24[3] = cheb12[3]+alam; 130 cheb24[23] = cheb12[3]-alam; 131 alam = x[6]*(v[2]-v[4]-v[6]); 132 cheb24[7] = cheb12[7]+alam; 133 cheb24[19] = cheb12[7]-alam; 134 alam = x[10]*v[2]-x[6]*v[4]+x[2]*v[6]; 135 cheb24[11] = cheb12[11]+alam; 136 cheb24[15] = cheb12[11]-alam; 137 for (i=1; i<=3; i++) { // end: 30 138 j = 8-i; 139 v[i] = fval[i]-fval[j]; 140 fval[i] = fval[i]+fval[j]; 141 l30: ;} 142 cheb12[5] = v[1]+x[8]*v[3]; 143 cheb12[9] = fval[1]-x[8]*fval[3]; 144 alam = x[4]*v[2]; 145 cheb24[5] = cheb12[5]+alam; 146 cheb24[21] = cheb12[5]-alam; 147 alam = x[8]*fval[2]-fval[4]; 148 cheb24[9] = cheb12[9]+alam; 149 cheb24[17] = cheb12[9]-alam; 150 cheb12[1] = fval[1]+fval[3]; 151 alam = fval[2]+fval[4]; 152 cheb24[1] = cheb12[1]+alam; 153 cheb24[25] = cheb12[1]-alam; 154 cheb12[13] = v[1]-v[3]; 155 cheb24[13] = cheb12[13]; 156 alam = 0.1e1/0.6e1; 157 for(i=2; i<=12; i++) { // end: 40 158 cheb12[i] = cheb12[i]*alam; 159 l40: ;} 160 alam = 0.5*alam; 161 cheb12[1] = cheb12[1]*alam; 162 cheb12[13] = cheb12[13]*alam; 163 for (i=2; i<=24; i++) { // end: 50 164 cheb24[i] = cheb24[i]*alam; 165 l50: ;} 166 cheb24[1] = 0.5*alam*cheb24[1]; 167 cheb24[25] = 0.5*alam*cheb24[25]; 168 return; 169 } 170 171 172 unittest 173 { 174 alias qcheb!float fqcheb; 175 alias qcheb!double dqcheb; 176 alias qcheb!real rqcheb; 177 }