/* * COPYRIGHT * * Copyright (C) 2005 Exstrom Laboratories LLC * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * A copy of the GNU General Public License is available on the internet at: * * http://www.gnu.org/copyleft/gpl.html * * or you can write to: * * The Free Software Foundation, Inc. * 675 Mass Ave * Cambridge, MA 02139, USA * * You can contact Exstrom Laboratories LLC via Email at: * * info(AT)exstrom.com * * or you can write to: * * Exstrom Laboratories LLC * P.O. Box 7651 * Longmont, CO 80501, USA * */ /* * gnm - generates the g(n,m) lattice Green function term using * equation (27) in cond-mat/0509002, and outputs the result * both symbolically and as a double. * To compile this program, you need to have a C compiler, the GMP (GNU MP Bignum) library, and the MPFR (multiple-precision floating-point computations with exact rounding) library. The GMP library can be gotten at http://www.swox.com/gmp/ and the MPFR library is at http://www.mpfr.org/ . The compiler that we use is gcc. Other compilers may or may not work. The program is compiled with the following command line: gcc -o gnm gnm.c -I/usr/local/include -L/usr/local/lib -lmpfr -lgmp Run the program with no parameters for details on using it. */ /* This program was adapted from the old program rnmd3. */ #include #include #include #include #include #include void s( signed long int n, signed long int m, mpz_t sz ); double gnm( signed long int n, signed long int m ); /***************************************************************************/ /* Implements the function: s(n,m) = p sum Hj j=0 where the Hj's are gotten from the recursion relation: Hj-1 = Hj * 4j*(n-j+1)*(2(n-j+1)+m)*(2j+m) ------------------------------ (n-2j+1)^2 * (n-2j+2)^2 and the initial value Hp is: Hp = { 1 n=even { (n+1)(n+m+1) n=odd Usage: s( signed long int n, signed long int m, mpz_t sz ); n = first parameter of function s(n,m) defined above. m = second parameter of function s(n,m) defined above. sz = the result of the function s(n,m) defined above. Note that if you want to change the value of a gmp data type within a function, you don't have to pass the address of the variable, since these variables are arrays of size 1, so just passing the variable and not its address is a call by reference. */ void s( signed long int n, signed long int m, mpz_t sz ) { signed long int i, p, q; mpz_t f; /* f is a dummy variable */ if( n < 0 ) { mpz_set_ui( sz, (unsigned long int)0 ); return; } /* initialize (allocate) dummy variable and set its value to zero */ mpz_init( f ); /* if n is odd, assign f=(n+1)*(n+m+1), otherwise f=1 */ if( n % 2 ) mpz_set_si(f,(n+1)*(n+m+1)); else mpz_set_si(f,(signed long int)1); /* set sz equal to f */ mpz_set( sz, f ); for( i = n/2; i > 0; --i ) { p = n-i+1; q = p-i; /* q = n-2*i+1 */ /* Multiply f by 4*i*p*(2*p+m)*(m+2*i). This is done in stages */ /* to avoid the multiplier being larger than a signed long int. */ mpz_mul_si( f, f, 4*i ); mpz_mul_si( f, f, p ); mpz_mul_si( f, f, 2*p+m ); mpz_mul_si( f, f, m+2*i ); /* Divide f by q*q*(q+1)*(q+1). This is done in stages to */ /* avoid the multiplier being larger than an unsigned long int. */ mpz_divexact_ui( f, f, (unsigned long int)q ); mpz_divexact_ui( f, f, (unsigned long int)q ); mpz_divexact_ui( f, f, (unsigned long int)(q+1) ); mpz_divexact_ui( f, f, (unsigned long int)(q+1) ); /* sz = sz + f */ mpz_add( sz, sz, f ); } /* clear dummy variable */ mpz_clear( f ); } /***************************************************************************/ /* Returns the lattice Green function value g(n,m) as a double, and also prints it out symbolically. */ double gnm( signed long int n, signed long int m ) { signed long int j, k; mpz_t d1, d2, d3; /* dn are dummy variables */ mpz_t sz; /* sz will hold the result of the s(n,k) function call */ mpq_t gnp; /* gnp = the rational number that is the non-pi part of the lattice Green function */ mpq_t gpi; /* gpi = the rational number that is the pi part of the lattice Green function */ mpq_t g; /* g = rational number used as a dummy variable */ mp_prec_t precision; /* precision = default precision of mpfr float types */ const unsigned long int sigfigswanted = 20; /* Number of sig figs wanted in final result (for double, 20 is good) */ size_t limbnumsize; /* the number of limbs [words] of the numerator of gnp */ unsigned long int numbits; /* the number of bits of the numerator of gnp */ mpfr_t gnp_mpf, gpi_mpf; /* the floating point version of the rationals Rnp and Rpi */ mpfr_t pi_mpf; /* pi_mpf = the transcendental number pi using the default precision. */ double gd; /* gd = gnp + gpi/pi converted to a double */ /* initialize gmp variables and set their values to zero */ mpz_init( sz ); mpz_init( d1 ); mpz_init( d2 ); mpz_init( d3 ); mpq_init( gnp ); mpq_init( gpi ); mpq_init( g ); /* k = integer part of (n+m)/2 */ k = ( n+m )/2; /* printf("Calculating the pi-part of the lattice Green function...\n"); */ for( j = 1; j <= k; ++j ) { /* This for loop calculates the sum in the pi part of the lattice Green function. */ /* get the value of the first s-function in the sum */ s( n+m-2*j, 2*j-2*m-1, sz ); /* set d1 = sz */ mpz_set( d1, sz ); /* If m != 0 then calculate the value of the second s-function */ /* in the sum, but if m = 0 then the 2nd s-function equals the */ /* 1st s-function, so don't recalculate it. */ if( m != 0 ) s( n-m-2*j, 2*j+2*m-1, sz ); /* set d2 = sz */ mpz_set( d2, sz ); /* d3 = d1 + d2 */ mpz_add( d3, d1, d2 ); /* if j is odd, multiply d3 by -1 */ if( (j % 2) == 1 ) mpz_neg( d3, d3 ); /* set the rational number g = d3 */ mpq_set_z( g, d3 ); /* divide g by 2*j-1 */ mpz_set_si( mpq_denref(g), 2*j-1 ); /* remove common factors in numerator and denominator of g */ mpq_canonicalize( g ); /* add g to gpi */ mpq_add( gpi, g, gpi ); /* remove common factors in numerator and denominator of gpi */ mpq_canonicalize( gpi ); } /* We still need to do a little more to finish calculating gpi */ /* if m is odd, multiply gpi by -1 */ if( (m % 2) == 1 ) mpq_neg( gpi, gpi ); /* gpi has now been calculated. Now moving on to gnp. */ /* printf("The pi-part has been calculated. Now moving on to the non-pi part...\n"); */ /* get the value of the s-function for gnp and store it in sz */ s( n-m-1, 2*m, sz ); /* set the rational number gnp = sz = integer */ mpq_set_z( gnp, sz ); /* divide gnp by 4 */ mpz_set_ui( mpq_denref(gnp), (unsigned long int)4 ); /* remove common factors in numerator and denominator of gnp */ mpq_canonicalize( gnp ); /* if m is odd, multiply gnp by -1 */ if( (m % 2) == 1 ) mpq_neg( gnp, gnp ); /* gnp has now been calculated. Lets print out the results. */ /* printf("The non-pi part has been calculated. The results are:\n"); */ /* print g = gnp + gpi */ printf("g(%d,%d) = ", n, m); mpq_out_str( stdout, 10, gnp ); printf(" + "); mpq_out_str( stdout, 10, gpi ); printf("/PI\n"); /* Now begin converting gnp + gpi/pi into a double here */ /* printf("Now beginning the conversion of gnp + gpi/pi into floating point...\n"); */ /* Set the default precision of mpfr_t type variables. */ /* We'll set the precision equal to the number of digits */ /* in the numerator of the rational form of Rnp, plus the */ /* number of sig figs (whole + fractional parts) wanted in */ /* the final result. */ /* First, get the size of the numerator of gnp in units of */ /* limbs, where a limb is a machine word. */ limbnumsize = mpz_size( mpq_numref(gnp) ); /* Convert the # of limbs of the gnp numerator to number of bits */ numbits = (unsigned long int) ( mp_bits_per_limb * ((signed long int)limbnumsize) ); /* Now we can set the default precision of mpfr_t type variables. */ precision = (mp_prec_t)numbits; /* Now add the number of sig figs wanted to the precision */ numbits = sigfigswanted * 53 / 15; /* A double has 53 bits per 15 sigfigs */ precision += numbits; /* precision = precision + (# of sig figs wanted in units of bits) */ mpfr_set_default_prec( precision ); /* Initialize the mpfr types gnp_mpf and gpi_mpf, and convert */ /* the rational type variables gnp & gpi into mpfr types. */ mpfr_init_set_q( gnp_mpf, gnp, GMP_RNDN ); mpfr_init_set_q( gpi_mpf, gpi, GMP_RNDN ); /* Set the variable pi_mpf equal to pi using the default precision. */ mpfr_init( pi_mpf ); mpfr_const_pi( pi_mpf, GMP_RNDN ); /* divide gpi_mpf by pi_mpf */ mpfr_div( gpi_mpf, gpi_mpf, pi_mpf, GMP_RNDN ); /* add gnp_mpf + gpi_mpf and store the result in gnp_mpf */ mpfr_add( gnp_mpf, gnp_mpf, gpi_mpf, GMP_RNDN ); /* Convert gnp_mpf, which now contains g(n,m), to a double */ gd = mpfr_get_d( gnp_mpf, GMP_RNDN ); /* printf("Finished the conversion of gnp + gpi/pi into floating point.\n"); */ /* free all the gmp variables */ mpz_clear( sz ); mpz_clear( d1 ); mpz_clear( d2 ); mpz_clear( d3 ); mpq_clear( gpi ); mpq_clear( gnp ); mpq_clear( g ); /* free all the mpfr variables */ mpfr_clear( gnp_mpf ); mpfr_clear( gpi_mpf ); mpfr_clear( pi_mpf ); return( gd ); } /***************************************************************************/ int main( int argc, char *argv[] ) { /* n and m are made signed integers to make operations with signed integers convenient */ signed long int n; /* n = x-axis coordinate of the lattice Green function, n >= 0. */ signed long int m; /* m = y-axis coordinate of the lattice Green function, m >= 0. */ double gd; /* gd = gnp + gpi/pi converted to a double */ if( argc < 3 ) { printf("\n gnm generates the g(n,m) lattice Green function term using"); printf("\n equation (27) in cond-mat/0509002, and outputs the result"); printf("\n both symbolically and as a double."); printf("\n Usage: gnm n m\n"); printf(" n = x-axis coordinate of the lattice Green function, n >= 0.\n"); printf(" m = y-axis coordinate of the lattice Green function, m >= 0.\n"); printf(" where n >= m.\n"); return(-1); } n = (signed long int)atoi(argv[1]); m = (signed long int)atoi(argv[2]); if( n < m ) { printf("ERROR: n must be greater than or equal to m.\n"); return(-1); } gd = gnm( n, m ); printf("g(%d,%d) = %1.15lf\n", n, m, gd); }