Can somebody suggest a code for calculating the nu function, as defined at http://mathworld.wolfram.com/NuFunction.html I am interested just in nu(x), not nu(x, alpha). Igor Volobouev ==== Possibly the following link to the following Matlab toolbox (free software under Gnu Public License) is of interest to you: http://ftp.cwi.nl/pauldz/Codes/LISQ/ It includes functions for the computation of moments and invariants (Hu), 2nd generation wavelet decomposition & reconstruction tools for images, visualization tools. Examples and documentation are included, e.g. see: http://ftp.cwi.nl/CWIreports/PNA/PNA-R0224.pdf The software is free under the terms of the GNU General Public License. Paul de Zeeuw Dr. Paul M. de Zeeuw CWI - PNA4 Amsterdam http://homepages.cwi.nl/~pauldz/ ==== In a book I have found the following formulation n xxxxxxxxxxxxxxxxxx x x x x x x x x x x xxxxxxxxxxxxxxxxxx j=1 It's look like a capital PI ! Can someone tell me waht its means ? -- Bernard Bour.8ee bernard@bouree.net ==== > > In a book I have found the following formulation > > > n > xxxxxxxxxxxxxxxxxx > x x > x x > x x > x x > x x > xxxxxxxxxxxxxxxxxx > j=1 > > It's look like a capital PI ! > > Can someone tell me waht its means ? > > > A capital 'pi' means product, just like a capital 'sigma' means summation. OUP ==== Dear all, does anyone know where one can find some c-source code performing a singular value decomposition of a complex matrix. Peter ==== > Dear all, > > does anyone know where one can find some c-source code performing a singular > value decomposition of a complex matrix. Lapack. Is written in fortran, but since you're only linking to it that doesn't matter. V. -- homepage: cs utk edu tilde lastname ==== If I remember correctly, the Fourier Transform (FT) of a convolution of two functions is equal to the product of the FT of the functions. If this is right, then I don't know why I'm having the following problem doing some experiments in Matlab. I'm convoluting two Gaussians, and I should get another Gaussian, but I don't: t=linspace(-10,10,1000); x=exp(-(t.^2)/2)/sqrt(2*pi); % Gaussian centered at 0 y=exp(-((t-1).^2)/2)/sqrt(2*pi); % Gaussian centered at 1 fx = fft(x); fy = fft(y); %corresponding FT's fconv = fx.*fy % product of FT's conv = ifft(fconv) % do an inverse FT However, when plotted, conv doesn't look like a Gaussian! Am I doing something wrong? Maybe I'm misunderstanding FFT's. Fernando G. del Cueto. ==== I am implementing a fast convolution program using FFTW library. I understand the constraints of the Discrete Convolution Theorem, and the issues of wrap around error and periodicity. Hence, I am taking the following approach at the moment. /*-------------------*/ Image: image[sizeX][sizeY] (Symmetrical)Filter: kernel[2k + 1][2k + 1] Step(1) Define padX = sizeX + k; padY = sizeY + k; Step (2) Pad image on two sides with zeros so that its size is padX * padY. I pad zeros on the right side and the bottom side of the image buffer. Step (3) Store kernel in a wrap around fashion with the size padX * padY. So, if Symmetrical (eg. Gaussian) kernel = [a, b, c d, e, f g, h, i] and padX = 8 and padY = 8 d, e, 0, 0, 0, 0, 0, f 0, 0, 0, 0, 0, 0, 0, 0 .... .... g,h, 0, 0, 0, 0, 0, i] Step (5) Do pixel by pixel complex multiplication of FFT output. Step (6) Take inverse FFT of the output Step (7) Read in the output Step (8) Remove zero paddings. /*---------*/ Well, this only works for small kernel sizes, because when the kernel size is bigger I get a shift in the blurred image. I tried various other ways but I haven't been successful in resolving this error. So, I ask for your help. The C++ code for the FFT functions is as follows: (1) /*-------Building Gaussian FFT--------------------*/ { float d; int counter; float normalise; int i; counter = 0; int y,x; normalise = 0; /////////////////////////////// //-2 -1 0 1 2 //-2 x //-1 x // 0 x // 1 x // 2 x /////////////////////////////// //e ^ -(x*x + y*y)/2 * sigma * sigma //////////////////////////////// //We only consider cases with k.kernelSide as odd for(y = -((k.kernelSide - 1)/2); y <= (k.kernelSide - 1)/2; y++) { for(x = -((k.kernelSide - 1)/2); x <= (k.kernelSide - 1)/2; x++) { d = (float)(x*x + y*y)/(float)(2.0f * k.sigma * k.sigma); kernel[counter] = pow(e,-d); normalise += kernel[counter]; counter++; } } for(i = 0; i < k.kernelSide*k.kernelSide; i++) kernel[i] = kernel[i]/normalise; } { float *kernelRow; float *kernelColumn; int y,x; kernelRow = new float[padX * k.kernelSide]; kernelColumn = new float[padSize]; /*---------Perform padding with zeros in wrap around fashion-------*/ kernel_FFT = new zomplex*[padY]; for(y = 0; y < padY; y++) { kernel_FFT[y] = new zomplex[padX]; for(x = 0; x < padX; x++) { kernel_FFT[y][x].im = 0; kernel_FFT[y][x].re = kernelColumn[x + y*padX]; } } ComputeFFT(kernel_FFT, -1); //Free the kernels that are no longer required delete [] kernelRow; delete [] kernelColumn; } /*--------------Wrap the kernel along rows-----------------*/ *kernelRow) { int y,x; for(y = 0; y < k.kernelSide; y++) { for(x = 0; x <= (k.kernelSide - 1)/2; x++) kernelRow[x + y*padX] = kernel[x + y*k.kernelSide]; for(x = (k.kernelSide - 1)/2 + 1; x < padX - (k.kernelSide - 1)/2; x++) kernelRow[x + y * padX] = 0; for(x = 1; x <= (k.kernelSide - 1)/2 ; x++) kernelRow[(padX - x) + y*padX] = kernel[(k.kernelSide - x) + y*k.kernelSide]; } } /*---------------Wrap the kernelRow along columns--------------*/ *kernelColumn) { int y,x; for(y = 0; y <= (k.kernelSide - 1)/2; y++) { for(x = 0; x < padX; x++) kernelColumn[x + y*padX] = kernel[x + y*padX]; } for(y = (k.kernelSide - 1)/2 + 1; y < padY - (k.kernelSide - 1)/2; y++) { for(x = 0; x < padX; x++) kernelColumn[x + y*padX] = 0; } for(y = 1; y <= (k.kernelSide - 1)/2; y++) { for(x = 0; x < padX; x++) kernelColumn[x + (padY - y) * padX] = kernel[x + (k.kernelSide - y)*padX]; } } (2) /*-----Function to compute FFT using FFTW library----------------*/ inline void ConvolutionOld::ComputeFFT(zomplex** array, int inv) { fftw_plan p; int y,x; fftw_complex *carray; carray = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * padSize); for(y = 0; y < padY; y++) { for(x = 0; x < padX; x++) { carray[x + y * padX][0] = array[y][x].re; carray[x + y * padX][1] = array[y][x].im; } } p = fftw_plan_dft_2d(padY, padX, carray, carray, inv, FFTW_ESTIMATE); fftw_execute(p); for(y = 0; y < padY; y++) { for(x = 0; x < padX; x++) { array[y][x].re = carray[x + y * padX][0]; array[y][x].im = carray[x + y * padX][1]; } } fftw_destroy_plan(p); fftw_free(carray); } (3)/*----------Function to compute image FFT-------------*/ /*----------Pad the image with zeros of the right and bottom--------------*/ inline void ConvolutionOld::ZeroPadImage() { int y,x; image_FFT = new zomplex*[padY]; for(y = 0; y < sizeY; y++) { image_FFT[y] = new zomplex[padX]; for(x = 0; x < sizeX; x++) { image_FFT[y][x].re = image[x + y*sizeX]; image_FFT[y][x].im = 0; } for(x = sizeX; x < padX; x++) { image_FFT[y][x].re = 0.0f; image_FFT[y][x].im = 0.0f; } } for(y = sizeY; y < padY; y++) { image_FFT[y] = new zomplex[padX]; for(x = 0; x < padX; x++) { image_FFT[y][x].re = 0.0f; image_FFT[y][x].im = 0.0f; } } ComputeFFT(image_FFT, -1); } (4) /*------Function to perform convolution and take inverse FFT----*/ inline void ConvolutionOld::Convolve(float *output) { int i,y,x; convolution_FFT = new zomplex*[padY]; //Do per pixel multiplication of the gaussian kernel and the image matrix for(y = 0; y < padY; y++) { convolution_FFT[y] = new zomplex[padX]; for(x = 0; x < padX; x++) { convolution_FFT[y][x].re = (image_FFT[y][x].re * kernel_FFT[y][x].re) - (image_FFT[y][x].im * kernel_FFT[y][x].im); convolution_FFT[y][x].im = (image_FFT[y][x].re * kernel_FFT[y][x].im) + (image_FFT[y][x].im * kernel_FFT[y][x].re); } } ComputeFFT(convolution_FFT, 1); //Find inverse FFT for(y = 0; y < sizeY; y++) { for(x = 0; x < sizeX; x++) output[x + y*sizeX] = convolution_FFT[y][x].re; } } This is a lot of code to post, but I am unable to understand the cause of the error. If someone has implemented 2D convolution using FFTW library and knows about how or how not do zero paddings, please help me. -Swati ==== could someone show me a way of obtaining the probability density function (pdf) of Z, where Z=arctan(Y/X), 0 < x < b, -b < y <0, and X and Y are independent, identically distributed. Also, if anybody knows of a good probability and random process book good for electrical engineers (not TOO much hard core math) where they have lots of examples that show more than just simple two-r.v. functions such as Z=X+Y, Z=X-Y, and Z=X/Y, I'm all ears. -Ed ==== why the area under a curve is given by its anti-derivative. Anybody should find this of interest. This link will take you to the spot in the appendix: www.precalculus.netfirms.com/#5 If someone you know is just learning the basics of calculus I could never retain the rules and methods of calculus without should help a new student make sense of it all. TR ==== There is a recent book (2002) on numerical methods with an accompanying numerical library that is little-known compared to Numerical Recipes but worth considering IMO. The book is by H. M. Antia, an astrophysicist from TIFR in Bombay, India, and it is published by Birkhauser. The Birkhauser and (more informative) Indian web sites are http://www.birkhauser.com/detail.tpl?cart=1074271936836359&isbn=3764367156 http://www.hindbook.com/own/antia.htm The book is hardbound, has 842 pages, and comes with a CD-ROM containing the Fortran 77 and C computer code and its documentation. The list price is $55. The book itself focuses more on numerical analysis and algorithms. The contents are as follows: 1. Introduction 2. Roundoff Error 3. Linear Algebraic Equations 4. Interpolation 5. Differentiation 6. Integration 7. Nonlinear Algebraic Equations 8. Optimisation 9. Functional Approximations 10. Algebraic Eigenvalue Problem 11. Ordinary Differential Equations 12. Integral Equations 13. Partial Differential Equations. Appendix B: Fortran Programs (On CD). Appendix C: C Programs (On CD). There are 223 Fortran subroutines and functions for chapters 2-13, each with a driver program and sample output. Each code is given in single, double, and quadruple-precision. Both Compaq Visual Fortran and Lahey/Fujitsu Fortran 95 can compile the single and double-precision versions, but only LF95 compiles the quadruple precision code. CVF 6.6 does not like (KIND=16), giving error messages like zroot.f(32) : Error: REAL(KIND=16) is not supported on this platform. [16] IMPLICIT REAL*16(A,B,D-H,O-Z) (Is there a work-around?). There are 210 C codes. Unlike Numerical Recipes, the source code in machine format is documented well. Also unlike NR, I did not see restrictions on distributing the source code. ==== > The book is hardbound, has 842 pages, and comes with a CD-ROM > containing the Fortran 77 and C computer code and its documentation. . . . . >Both Compaq Visual Fortran > and Lahey/Fujitsu Fortran 95 can compile the single and > double-precision versions, but only LF95 compiles the quadruple > precision code. CVF 6.6 does not like (KIND=16), giving error messages Doesn't sound like Fortran 77 to me. But can you characterize the programs in more detail? What can you find in this book that is not adequately treated in Numerical Recipes? ==== >> Is there a relationship known for the distance between the largest root >> of an orthogonal polynomial and the endpoint of the orthogonality >> interval? I would like to know how fast this distance shrinks as a >> function of the degree of the polynomial... > > > > 1) Chebyshev 1st kind, orthogonal on [-1,1]: x_n = cos((2n-1)pi/2n) > > 2) Legendre, orthogonal on [-1,1]: x_n = 1-(j_n)^2/2 * 1/n^2 + O(1/n^3). > Here j_n is the nth biggest positive zero of Bessel function J_0(x). > > All the best, > Teijo > > non-classical orthogonal polynomials at hand. I have the orthogonality > interval and the weight function, but not a closed form of the polynomials. > Is there some general theory? I've just ordered Szeg.9a's book Orthogonal > Polynomials. I hope to find something in there.... > > Bye, > gert the zeros of orthogonal polynomial are analyzed in detail in the mentioned book from Szeg.9a for ultraspherical Jacobi polynomials (orthogonal wrt. the weight function rho(x)=(1-x)^alpha (1+x)^beta on [-1,1], -0.5 You know how in base 10, if the sum of the digits of any number add up > to a multiple of 3 or 9, then that number is not prime? ... I've perhaps slipped a cog, but the sum of the digits of, say, 18 add up to a multiple of both 3 and 9, and yet 18 isn't prime. What have I missed? -- Alan line, otherwise I won't see it. Ain't spam a wonderful thing? ==== David Blume: > You know how in base 10, if the sum of the digits of any number add up > to a multiple of 3 or 9, then that number is not prime? ... Counterexample: 3 is prime. But we knew what you meant. Alan Waldock: > I've perhaps slipped a cog, but the sum of the digits of, say, 18 add up to > a multiple of both 3 and 9, and yet 18 isn't prime. What have I missed? The second-last word of David's posting? -- Mark Brader, Toronto, msb@vex.net Irving Thalberg's advice on GONE WITH THE WIND: Forget it, Louis. No Civil War picture ever made a nickel. ==== > David Blume: > You know how in base 10, if the sum of the digits of any number add up > to a multiple of 3 or 9, then that number is not prime? ... Counterexample: 3 is prime. But we knew what you meant. Alan Waldock: > I've perhaps slipped a cog, but the sum of the digits of, say, 18 add up to > a multiple of both 3 and 9, and yet 18 isn't prime. What have I missed? The second-last word of David's posting? Colour me foolish. Mustn't post so late at night!! -- Alan line, otherwise I won't see it. Ain't spam a wonderful thing? ==== > You know how in base 10, if the sum of the digits of any number add up > to a multiple of 3 or 9, then that number is not prime? Can it be > proven that it works in the general case? I know it to be true, but > don't know the proof. > > That is, in any base b, if the sum of the digits of any positive > number n add up to a multiple of any of the factors of (b - 1), then > that number is not prime. > > For example, in base 241, there is no prime number whose sum of the > digits add up to multiples of the digits represented by 2, 3, 4, 5, 6, > 8, 10, 12, 15, 16, 20, 24, 30, 40, 48, 60, 80, or 240. (These digits > were written in base-10 for simplicity's sake. But they are indeed > single digits in base 241.) > > Ex., In base 241, the prime 65393 would be written 1,30,82. (Again, > digits in base 10 for simplicity.) The sum of those digits, 113, > isn't divisible by any of the factors of 240. If n divides b-1, then b = 1 mod n, hence by an easy induction b^k = 1 mod n, (k >= 0) hence sum_k c_k b^k = sum_k c_k mod n. You do have to be careful about niggling cases, e.g. b = 3 and n = 2. See if YOU can find the most general niggling case. --Ron Bruck ==== ) You know how in base 10, if the sum of the digits of any number add up ) to a multiple of 3 or 9, then that number is not prime? Can it be ) proven that it works in the general case? I know it to be true, but ) don't know the proof. ) ) That is, in any base b, if the sum of the digits of any positive ) number n add up to a multiple of any of the factors of (b - 1), then ) that number is not prime. ) ) For example, in base 241, there is no prime number whose sum of the ) digits add up to multiples of the digits represented by 2, 3, 4, 5, 6, ) 8, 10, 12, 15, 16, 20, 24, 30, 40, 48, 60, 80, or 240. (These digits ) were written in base-10 for simplicity's sake. But they are indeed ) single digits in base 241.) ) ) Ex., In base 241, the prime 65393 would be written 1,30,82. (Again, ) digits in base 10 for simplicity.) The sum of those digits, 113, ) isn't divisible by any of the factors of 240. Yeah, that's easy. A number in base B is equal to the sum of the digits, where each digit is multiplied by an appropriate power of B. If F is a factor of B-1, then the following equations hold: X * (B-1) = 0, modulo F. X * B = (X * (B-1)) + X = 0 + X = X, modulo F. X*B^N = X*B^N-1, modulo F. X*B^N = X, modulo F. D1*B^0 + D2*B^1 + D3*B^2 +...+ Dn * B^n = D1 + D2 + D3 +...+ Dn, modulo F. IOW: the sum of the digits of a number is equal to that number, modulo F. A number is only divisible by F if it is equal to 0, modulo F. SaSW, Willem -- made in the above text. For all I know I might be drugged or something.. No I'm not paranoid. You all think I'm paranoid, don't you ! #EOT X-Received: (from approve@localhost) by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision: 1.9 primary) id i0EESHc04831; Wed, 14 Jan 2004 09:28:17 -0500 ==== on;0208-541-1735....thanks,kevin. X-Received: (from approve@localhost) by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision: 1.9 primary) id i0EJROH28736; ==== I am visual.net and I successfully compiled the lapack library. However when I introduce it to my code I have some library conflict such as: Linking... MSVCRTD.lib(MSVCR71D.dll) : error LNK2005: _printf already defined in LIBCD.lib(printf.obj) MSVCRTD.lib(MSVCR71D.dll) : error LNK2005: _free already defined in LIBCD.lib(dbgheap.obj) MSVCRTD.lib(MSVCR71D.dll) : error LNK2005: _fprintf already defined in LIBCD.lib(fprintf.obj) MSVCRTD.lib(MSVCR71D.dll) : error LNK2005: _malloc already defined in LIBCD.lib(dbgheap.obj) MSVCRTD.lib(MSVCR71D.dll) : error LNK2005: _exit already defined in LIBCD.lib(crt0dat.obj) MSVCRTD.lib(ti_inst.obj) : error LNK2005: private: __thiscall type_info::type_info(class type_info const &) (??0type_info@@AAE@ABV0@@Z) already defined in LIBCD.lib(typinfo.obj) MSVCRTD.lib(ti_inst.obj) : error LNK2005: private: class type_info & __thiscall type_info::operator=(class type_info const &) (??4type_info@@AAEAAV0@ABV0@@Z) already defined in LIBCD.lib(typinfo.obj) LINK : warning LNK4098: defaultlib 'MSVCRTD' conflicts with use of other libs; use /NODEFAULTLIB:library Debug/PP_stat.exe : fatal error LNK1169: one or more multiply defined symbols found Do you why is that? What should I do? David simonin X-Received: (from approve@localhost) by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision: 1.9 primary) id i0G0SWk29958; ==== I have indexed only the non-zero entries of a sparse matrix using the expression: ij=(row-1)*matrix_size +columns. Is it possible to reconstruct the the row and the columns of the matrix depending only the index, ij. Where matrix_size is the size of a given square matrix. This is becuase I can not save the row and the columns of each entry in the first pass of my program. I have tried to find the row and the columns like this: columns=ij-matrix_size*INT(ij/matrix_size) row=INT((ij-columns)/matrix_size)+columns, where INT produce the integer part of the expressions above. I, however, did not quite get the right answer ? Is there any one out there who can find a solution for this ? Yacob