Changeset 5529:4f7db61caccf


Ignore:
Timestamp:
07/30/07 01:14:12 (6 years ago)
Author:
William Stein <wstein@…>
Branch:
default
Message:

New version of number of partitions code from J Bober.

Location:
sage/combinat
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • sage/combinat/combinat.py

    r5525 r5529  
    11771177    return eval(ans.replace('\n','')) 
    11781178 
    1179 def number_of_partitions(n,k=None, algorithm='gap'): 
     1179def number_of_partitions(n,k=None, algorithm='bober'): 
    11801180    r""" 
    11811181    Returns the size of partitions_list(n,k). 
    … …  
    11861186             cardinality of the set of all (unordered) partitions of 
    11871187             the positive integer n into sums with k summands. 
    1188         algorithm -- (default: 'gap') 
    1189             'gap' -- use GAP (VERY *slow*) 
     1188        algorithm -- (default: 'bober') 
    11901189            'bober' -- use Jonathon Bober's implementation (*very* fast, 
    11911190                      but new and not well tested yet). 
     1191            'gap' -- use GAP (VERY *slow*) 
    11921192            'pari' -- use PARI.  Speed seems the same as GAP until $n$ is 
    11931193                      in the thousands, in which case PARI is faster. *But* 
    … …  
    12111211        sage: len(v) 
    12121212        7 
    1213         sage: number_of_partitions(5) 
     1213        sage: number_of_partitions(5, algorithm='gap') 
    12141214        7 
    12151215        sage: number_of_partitions(5, algorithm='pari') 

  • sage/combinat/partitions_c.cc

    r5526 r5529  
    1 /* 
    2     Program to compute the number of partitions of n. 
    3      
    4     Author: Jonathan Bober 
    5     Version: .1 (7/28/2007) 
    6  
    7     This program is free software; you can redistribute it and/or modify 
    8     it under the terms of the GNU General Public License as published by 
    9     the Free Software Foundation; either version 2 of the License, or 
    10     (at your option) any later version. 
    11  
    12     This program is distributed in the hope that it will be useful, 
    13     but WITHOUT ANY WARRANTY; without even the implied warranty of 
    14     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
    15     GNU General Public License for more details. 
    16  
    17     You should have received a copy of the GNU General Public License 
    18     along with this program.  If not, see <http://www.gnu.org/licenses/>. 
    19 */ 
     1/*      Author:     Jonathan Bober 
     2 *      Version:    .3 
     3 * 
     4 *      This program computes p(n), the number of integer partitions of n, using Rademacher's 
     5 *      formula. (See Hans Rademacher, On the Partition Function p(n), 
     6 *      Proceedings of the London Mathematical Society 1938 s2-43(4):241-254; doi:10.1112/plms/s2-43.4.241, 
     7 *      currently at 
     8 * 
     9 *      http://plms.oxfordjournals.org/cgi/content/citation/s2-43/4/241 
     10 * 
     11 *      if you have access.) 
     12 * 
     13 *      We use the following notation: 
     14 * 
     15 *      p(n) = lim_{n --> oo} t(n,N) 
     16 * 
     17 *      where 
     18 * 
     19 *      t(n,N) = sum_{k=1}^N a(n,k) f_n(k), 
     20 * 
     21 *      where 
     22 * 
     23 *      a(n,k) = sum_{h=1, (h,k) = 1}^k exp(\pi i s(h,k) - 2 \pi i h  n / k) 
     24 * 
     25 *      and 
     26 *       
     27 *      f_n(k) = \pi sqrt{2} cosh(A_n/(sqrt{3}*k))/(B_n*k) - sinh(C_n/k)/D_n; 
     28 * 
     29 *      where 
     30 * 
     31 *      s(h,k) = \sum_{j=1}^{k-1}(j/k)((hj/k)) 
     32 * 
     33 *      A_n = sqrt{2} \pi * sqrt{n - 1/24} 
     34 *      B_n = 2 * sqrt{3} * (n - 1/24) 
     35 *      C_n = sqrt{2} * \pi * sqrt{n - 1.0/24.0} / sqrt{3} 
     36 *      D_n = 2 * (n - 1/24) * sqrt{n - 1.0/24.0} 
     37 * 
     38 *      and, finally, where ((x)) is the sawtooth function ((x)) = {x} - 1/2 if x is not an integer, 0 otherwise. 
     39 * 
     40 *      Some clever tricks are used in the computation of s(h,k), and perhaps at least 
     41 *      some of these are due to Apostol. (I don't know a reference for this.) 
     42 * 
     43 *      TODO: fix memory leaks 
     44 *          -- All mpfr_t, mpq_t, and mpz_t variables should be cleared before the function 
     45 *             mp_t(n) returns, but currently this is not done. This doesn't really matter 
     46 *             when this is a standalone program, but will matter if this is called as 
     47 *             a library function from another program. 
     48 * 
     49 *          -- Search source code for other TODO comments. 
     50 * 
     51 *      OTHER CREDITS: 
     52 *       
     53 *      I looked source code written by Ralf Stephan, currently available at 
     54 * 
     55 *              http://www.ark.in-berlin.de/part.c 
     56 * 
     57 *      while writing this code, but didn't really make use of it, except for the 
     58 *      function cospi(), currently included in the source (slightly modified), but not 
     59 *      used. 
     60 * 
     61 *      More useful were notes currently available at 
     62 * 
     63 *              http://www.ark.in-berlin.de/part.pdf 
     64 * 
     65 *      and others at 
     66 * 
     67 *              http://www.math.uwaterloo.ca/~dmjackso/CO630/ptnform1.pdf 
     68 * 
     69 *      Also, it is worth noting that the code for GCD(), while trivial, 
     70 *      was directly copied from the NTL source code. 
     71 * 
     72 *      Also, Bill Hart made some comments about ways to speed up this computation on the SAGE 
     73 *      mailing list. 
     74 * 
     75 *      This program is free software; you can redistribute it and/or modify 
     76 *      it under the terms of the GNU General Public License as published by 
     77 *      the Free Software Foundation; either version 2 of the License, or 
     78 *      (at your option) any later version. 
     79 * 
     80 *      This program is distributed in the hope that it will be useful, 
     81 *      but WITHOUT ANY WARRANTY; without even the implied warranty of 
     82 *      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
     83 *      GNU General Public License for more details. 
     84 * 
     85 *      You should have received a copy of the GNU General Public License 
     86 *      along with this program.  If not, see <http://www.gnu.org/licenses/>.  
     87 */ 
     88 
     89 
     90 
     91 
     92 
    2093#include <stdio.h> 
    2194 
    … …  
    33106 
    34107const bool debugs = false; 
    35  
    36108//const bool debugf = true; 
    37109const bool debugf = false; 
    … …  
    47119const unsigned int min_precision = 53; 
    48120 
     121bool test(bool longtest = false); 
     122 
     123void cospi (mpfr_t res, mpfr_t x); 
     124 
     125 
     126// The following function can be useful for debugging in come circumstances, but should not be used for anything else  
     127// unless it is rewritten. 
     128int grab_last_digits(char * output, int n, mpfr_t x) { 
     129    // fill output with the n digits of x that occur 
     130    // just before the decimal point 
     131    // Note: this assumes that x has enough digits and enough 
     132    // precision -- otherwise bad things can happen 
     133    // 
     134    // returns: the number of digits to the right of the decimal point 
     135 
     136    char * temp; 
     137    mp_exp_t e; 
     138 
     139    temp = mpfr_get_str(NULL, &e, 10, 0, x, GMP_RNDN); 
     140 
     141    int retval; 
     142 
     143    if(e > 0) { 
     144        strncpy(output, temp + e - n, n);     
     145        retval =  strlen(temp + e); 
     146    } 
     147    else { 
     148        for(int i = 0; i < n; i++) 
     149            output[i] = '0'; 
     150        retval = strlen(temp); 
     151    } 
     152    output[n] = '\0'; 
     153 
     154    mpfr_free_str(temp); 
     155 
     156    return retval; 
     157} 
     158 
     159 
    49160long GCD(long a, long b); 
    50161 
    … …  
    54165 
    55166mpz_t ztemp1, ztemp2, ztemp3; 
    56  
    57 mpfr_t mp_one_over_12, mp_one_over_24, mp_sqrt2, mp_sqrt3, mp_pi; 
     167mpq_t qtemp1, qtemp2, qtemp3; 
     168 
     169mpfr_t mp_one_over_12, mp_one_over_24, mp_sqrt2, mp_sqrt3, mp_pi, half, fourth; 
    58170mpfr_t mp_A, mp_B, mp_C, mp_D; 
    59171mp_rnd_t round_mode = GMP_RNDN; 
    60172 
    61173mpfr_t tempa1, tempa2, tempf1, tempf2, temps1, temps2, tempt1, tempt2;  // temp variables for different functions, with precision set and cleared by mp_set_precision 
     174mpfr_t tempc1, tempc2; // temp variable used by cospi() 
     175 
    62176 
    63177bool mp_vars_initialized = false; 
    … …  
    65179mp_prec_t mp_precision; 
    66180 
    67 void mp_init(unsigned int prec, unsigned int n); 
     181void initialize_constants(unsigned int prec, unsigned int n); 
    68182//void mp_f_precompute(unsigned int n); 
    69183void mp_f(mpfr_t result, unsigned int k); 
    … …  
    85199//low precision (double) versions of the functions: 
    86200 
    87 double df_A, df_B, df_C, df_D; 
    88  
    89 void d_f_precompute(unsigned int n); 
     201double d_A, d_B, d_C, d_D; 
     202 
    90203double d_f(unsigned int k); 
    91204double d_s(unsigned int h,unsigned int q); 
    92 double d_A(unsigned int n, unsigned int k); 
     205double d_a(unsigned int n, unsigned int k); 
    93206double d_t(unsigned int n, unsigned int N); 
    94207 
    … …  
    115228    // to be sure that we compute the correct answer 
    116229 
    117     /* I don't think that we need to be that careful about this. 
    118      * 
    119     mpfr_t bits; 
    120     mpfr_t t1;                                          // we refer to t1 as temp1 
    121      
    122     mpfr_init2(bits, 32);     // use 32 bits of precision to compute the approximation 
    123  
    124     mpfr_init2(t1, 32); 
    125  
    126     //    Error = exp( pi* sqrt(2n/3)); 
    127  
    128     mpfr_set_ui(t1, 2, round_mode);                     // temp1 = 2 
    129     mpfr_mul_ui(t1, t1, n, round_mode);                 // temp1 = 2n 
    130     mpfr_div_ui(t1, t1, 3, round_mode);                 // temp1 = 2n/3 
    131     mpfr_sqrt(t1, t1, round_mode);                      // temp1 = sqrt(2n/3) 
    132  
    133     mpfr_const_pi(bits, round_mode);                    // bits = pi 
    134  
    135     mpfr_mul(bits, bits, t1, round_mode);               // bits = pi * sqrt(2n/3) 
    136  
    137     mpfr_const_log2(t1, round_mode);                    // temp1 = log(2) 
    138     mpfr_div(bits, bits, t1, round_mode);               // bits = pi * sqrt(2n/3)/log(2) 
    139     */ 
    140230    unsigned int result = (unsigned int)(ceil(3.1415926535897931 * sqrt(2.0 * double(n)/ 3.0) / log(2))) + 3; 
    141     if(debug) 
    142         cout << "Using initial precision of " << result << " bits." << endl; 
     231    if(debug) cout << "Using initial precision of " << result << " bits." << endl; 
     232     
    143233    if(result > min_precision) { 
    144234        return result; 
    145235    } 
     236     
    146237    else return min_precision; 
    147238 
    … …  
    152243    //      log(A/sqrt(N) + B*sqrt(N/(n-1))*sinh(C * sqrt(n) / N) / log(2) 
    153244    // 
    154     //  maybe there is a better way... 
    155  
    156     // cout << N; 
     245    //  where A, B, and C are the constants listed below. These error bounds 
     246    //  are given in the paper by Rademacher listed at the top of this file. 
     247    // 
    157248 
    158249    mpfr_t A, B, C; 
    … …  
    189280    mpfr_add(error, error, t1, round_mode); // error = (ERROR ESTIMATE) 
    190281     
    191     unsigned int p = mpfr_get_exp(error) + 3;   // I am almost certain that this does the right thing 
     282    unsigned int p = mpfr_get_exp(error);   // I am almost certain that this does the right thing 
    192283                                                // (The 3 is for good luck.) 
    193284     
    … …  
    210301 
    211302 
    212  
    213  
    214303    if(p > min_precision) { 
    215         return p;                           // don't want to return < 32 
    216                                             // (for now, at least -- we can be more careful) 
     304        return p;                           // don't want to return < min_precision 
     305                                            // Note that when we hit the minimum precision 
     306                                            // we should switch over to using C doubles instead 
     307                                            // of mpfr types. 
    217308    } 
    218309    return min_precision; 
    … …  
    220311 
    221312double compute_remainder(unsigned int n, unsigned int N) { 
     313    // This computes the remainer left after N terms have been computed. 
     314    // The formula is exactly the same as the one used to compute the required 
     315    // precision, but once we know the necessary precision is small, we can 
     316    // call this function to determine the actual error (rather than the precision). 
     317    // 
     318    // Generally, this is only called once we know that the necessary 
     319    // precision is <= min_precision, because then the error is small 
     320    // enough to fit into a double, and also, we know that we are 
     321    // getting close to the correct answer. 
     322 
    222323    double A = 1.11431833485164; 
    223324    double B = 0.059238439175445; 
    … …  
    230331 
    231332void mp_set_precision(unsigned int prec) { 
     333    // 
     334    // Clear and initialize some "temp" variables that are used in the computation of various functions. 
     335    // 
     336    // We do this in this auxilliary function (and endure the pain of 
     337    // the extra global variables) so that we only have to initialize/clear 
     338    // these variables once every time the precision changes. 
     339    // 
     340    // NAMING CONVENTIONS: 
     341    // 
     342    // -tempa1 and tempa2 are two variables available for use 
     343    //      in the function mp_a() 
     344    // 
     345    // -tempf1 and tempf2 are two variables available for use 
     346    //      in the function mp_f() 
     347    // 
     348    // -etc... 
     349 
    232350    static bool init = false; 
    233351    if(init) { 
    … …  
    240358        mpfr_clear(temps1); 
    241359        mpfr_clear(temps2); 
     360 
     361        mpfr_clear(tempc1); 
     362        mpfr_clear(tempc2); 
    242363    } 
    243364    mpfr_init2(tempa1, prec); 
    … …  
    249370    mpfr_init2(temps1, prec); 
    250371    mpfr_init2(temps2, prec); 
     372    mpfr_init2(tempc1, prec); 
     373    mpfr_init2(tempc2, prec); 
     374 
    251375 
    252376    init = true; 
    253377} 
    254378 
    255 void mp_init(unsigned int prec, unsigned int n) { 
     379 
     380void initialize_constants(unsigned int prec, unsigned int n) { 
     381    // The variables mp_A, mp_B, mp_C, and mp_D are used for 
     382    // A_n, B_n, C_n, and D_n listed at the top of this file. 
     383    // 
     384    // They depend only on n, so we compute them just once in this function, 
     385    // and then use them many times elsewhere. 
     386    // 
     387    // Also, we precompute some extra constants that we use a lot, such as 
     388    // sqrt2, sqrt3, pi, 1/24, 1/12, etc. 
    256389    static bool init = false; 
    257390    mp_precision = prec; 
    … …  
    260393    if(init) { 
    261394        mpfr_clear(mp_one_over_12); mpfr_clear(mp_one_over_24); mpfr_clear(mp_sqrt2); mpfr_clear(mp_sqrt3); mpfr_clear(mp_pi); 
    262         mpfr_clear(mp_A); mpfr_clear(mp_B); mpfr_clear(mp_C); mpfr_clear(mp_D); 
     395        mpfr_clear(mp_A); mpfr_clear(mp_B); mpfr_clear(mp_C); mpfr_clear(mp_D); mpfr_clear(half); mpfr_clear(fourth); 
    263396 
    264397        mpz_clear(ztemp1); 
    265398        mpz_clear(ztemp2); 
    266399        mpz_clear(ztemp3); 
     400         
     401        mpq_clear(qtemp1); 
     402        mpq_clear(qtemp2); 
     403        mpq_clear(qtemp3); 
    267404    } 
    268405    mpfr_init2(mp_one_over_12,p); mpfr_init2(mp_one_over_24,p); mpfr_init2(mp_sqrt2,p); mpfr_init2(mp_sqrt3,p); mpfr_init2(mp_pi,p); 
    269     mpfr_init2(mp_A,p); mpfr_init2(mp_B,p); mpfr_init2(mp_C,p); mpfr_init2(mp_D,p); 
     406    mpfr_init2(mp_A,p); mpfr_init2(mp_B,p); mpfr_init2(mp_C,p); mpfr_init2(mp_D,p); mpfr_init2(fourth, p); mpfr_init2(half, p); 
    270407    mpz_init(ztemp1); 
    271408    mpz_init(ztemp2); 
    272409    mpz_init(ztemp3); 
    273410 
     411    mpq_init(qtemp1); 
     412    mpq_init(qtemp2); 
     413    mpq_init(qtemp3); 
     414     
    274415    init = true; 
    275416     
    276     mpfr_set_ui(mp_one_over_12, 1, round_mode);                         // mp_one_over_12 = 1/12 
    277     mpfr_div_ui(mp_one_over_12, mp_one_over_12, 12, round_mode);        // 
    278  
    279     mpfr_set_ui(mp_one_over_24, 1, round_mode);                         // mp_one_over_24 = 1/24 
    280     mpfr_div_ui(mp_one_over_24, mp_one_over_24, 24, round_mode);        // 
    281  
    282     mpfr_t n_minus;                                                     // 
    283     mpfr_init2(n_minus, p);                                             // 
    284     mpfr_set_ui(n_minus, n, round_mode);                                // n_minus = n 
    285     mpfr_sub(n_minus, n_minus, mp_one_over_24,round_mode);              // n_minus = n - 1/24 
    286  
    287     mpfr_t sqrt_n_minus;                                                // 
    288     mpfr_init2(sqrt_n_minus, p);                                        // 
    289     mpfr_sqrt(sqrt_n_minus, n_minus, round_mode);                       // n_minus = sqrt(n-1/24) 
    290  
    291  
    292     mpfr_sqrt_ui(mp_sqrt2, 2, round_mode);                              // mp_sqrt2 = sqrt(2) 
    293     mpfr_sqrt_ui(mp_sqrt3, 3, round_mode);                              // mp_sqrt3 = sqrt(3) 
    294     mpfr_const_pi(mp_pi, round_mode);                                   // mp_pi = pi 
    295  
    296     //mp_A = sqrt(2) * 3.1415926535897931 * sqrt(n - 1.0/24.0);----------- 
    297                                                                         // 
    298     mpfr_set(mp_A, mp_sqrt2, round_mode);                               // mp_A = sqrt(2) 
    299     mpfr_mul(mp_A, mp_A, mp_pi, round_mode);                            // mp_A = sqrt(2) * pi 
    300     mpfr_mul(mp_A, mp_A, sqrt_n_minus, round_mode);                     // mp_A = sqrt(2) * pi * sqrt(n - 1/24) 
    301     //-------------------------------------------------------------------- 
    302  
    303     //cout << "n_minus_1/24: "; 
    304     //mpfr_out_str(stdout, 10, 20, n_minus, round_mode); 
    305     //cout << endl; 
    306  
    307     //mp_B = 2.0 * sqrt(3) * (n - 1.0/24.0);------------------------------ 
    308     mpfr_set_ui(mp_B, 2, round_mode);                                   // mp_A = 2 
    309     mpfr_mul(mp_B, mp_B, mp_sqrt3, round_mode);                         // mp_A = 2*sqrt(3) 
    310     mpfr_mul(mp_B, mp_B, n_minus, round_mode);                          // mp_A = 2*sqrt(3)*(n-1/24) 
    311     //-------------------------------------------------------------------- 
    312  
    313     //mp_C = sqrt(2) * pi * sqrt(n - 1.0/24.0) / sqrt(3);----------------- 
    314     mpfr_set(mp_C, mp_sqrt2, round_mode);                               // mp_C = sqrt(2) 
    315     mpfr_mul(mp_C, mp_C, mp_pi, round_mode);                                  // mp_C = sqrt(2) * pi 
    316     mpfr_mul(mp_C, mp_C, sqrt_n_minus, round_mode);                           // mp_C = sqrt(2) * pi * sqrt(n - 1/24) 
    317     mpfr_div(mp_C, mp_C, mp_sqrt3, round_mode);                               // mp_C = sqrt(2) * pi * sqrt(n - 1/24) / sqrt3 
    318     //-------------------------------------------------------------------- 
    319  
    320     //mp_D = 2.0 * (n - 1.0/24.0) * sqrt(n - 1.0/24.0);------------------- 
    321     mpfr_set_ui(mp_D, 2, round_mode);                                   // mp_D = 2 
    322     mpfr_mul(mp_D, mp_D, n_minus, round_mode);                          // mp_D = 2 * (n - 1/24) 
    323     mpfr_mul(mp_D, mp_D, sqrt_n_minus, round_mode);                     // mp_D = 2 * (n - 1/24) * sqrt(n - 1/24) 
    324     //-------------------------------------------------------------------- 
    325      
     417    mpfr_set_ui(mp_one_over_12, 1, round_mode);                             // mp_one_over_12 = 1/12 
     418    mpfr_div_ui(mp_one_over_12, mp_one_over_12, 12, round_mode);            // 
     419 
     420    mpfr_set_ui(mp_one_over_24, 1, round_mode);                             // mp_one_over_24 = 1/24 
     421    mpfr_div_ui(mp_one_over_24, mp_one_over_24, 24, round_mode);            // 
     422 
     423    mpfr_set_ui(half, 1, round_mode);                                       // 
     424    mpfr_div_ui(half, half, 2, round_mode);                                    // half = 1/2 
     425    mpfr_div_ui(fourth, half, 2, round_mode);                                  // fourth = 1/4 
     426 
     427    mpfr_t n_minus;                                                         // 
     428    mpfr_init2(n_minus, p);                                                 // 
     429    mpfr_set_ui(n_minus, n, round_mode);                                    // n_minus = n 
     430    mpfr_sub(n_minus, n_minus, mp_one_over_24,round_mode);                  // n_minus = n - 1/24 
     431 
     432    mpfr_t sqrt_n_minus;                                                    // 
     433    mpfr_init2(sqrt_n_minus, p);                                            // 
     434    mpfr_sqrt(sqrt_n_minus, n_minus, round_mode);                           // n_minus = sqrt(n-1/24) 
     435 
     436 
     437    mpfr_sqrt_ui(mp_sqrt2, 2, round_mode);                                  // mp_sqrt2 = sqrt(2) 
     438    mpfr_sqrt_ui(mp_sqrt3, 3, round_mode);                                  // mp_sqrt3 = sqrt(3) 
     439    mpfr_const_pi(mp_pi, round_mode);                                       // mp_pi = pi 
     440 
     441    //mp_A = sqrt(2) * 3.1415926535897931 * sqrt(n - 1.0/24.0);--------------- 
     442                                                                            // 
     443    mpfr_set(mp_A, mp_sqrt2, round_mode);                                   // mp_A = sqrt(2) 
     444    mpfr_mul(mp_A, mp_A, mp_pi, round_mode);                                // mp_A = sqrt(2) * pi 
     445    mpfr_mul(mp_A, mp_A, sqrt_n_minus, round_mode);                         // mp_A = sqrt(2) * pi * sqrt(n - 1/24) 
     446    //------------------------------------------------------------------------ 
     447 
     448    //mp_B = 2.0 * sqrt(3) * (n - 1.0/24.0);---------------------------------- 
     449    mpfr_set_ui(mp_B, 2, round_mode);                                       // mp_A = 2 
     450    mpfr_mul(mp_B, mp_B, mp_sqrt3, round_mode);                             // mp_A = 2*sqrt(3) 
     451    mpfr_mul(mp_B, mp_B, n_minus, round_mode);                              // mp_A = 2*sqrt(3)*(n-1/24) 
     452    //------------------------------------------------------------------------ 
     453 
     454    //mp_C = sqrt(2) * pi * sqrt(n - 1.0/24.0) / sqrt(3);--------------------- 
     455    mpfr_set(mp_C, mp_sqrt2, round_mode);                                   // mp_C = sqrt(2) 
     456    mpfr_mul(mp_C, mp_C, mp_pi, round_mode);                                // mp_C = sqrt(2) * pi 
     457    mpfr_mul(mp_C, mp_C, sqrt_n_minus, round_mode);                         // mp_C = sqrt(2) * pi * sqrt(n - 1/24) 
     458    mpfr_div(mp_C, mp_C, mp_sqrt3, round_mode);                             // mp_C = sqrt(2) * pi * sqrt(n - 1/24) / sqrt3 
     459    //------------------------------------------------------------------------ 
     460 
     461    //mp_D = 2.0 * (n - 1.0/24.0) * sqrt(n - 1.0/24.0);----------------------- 
     462    mpfr_set_ui(mp_D, 2, round_mode);                                       // mp_D = 2 
     463    mpfr_mul(mp_D, mp_D, n_minus, round_mode);                              // mp_D = 2 * (n - 1/24) 
     464    mpfr_mul(mp_D, mp_D, sqrt_n_minus, round_mode);                         // mp_D = 2 * (n - 1/24) * sqrt(n - 1/24) 
     465    //------------------------------------------------------------------------ 
     466     
     467 
    326468    mpfr_clear(n_minus); 
    327469    mpfr_clear(sqrt_n_minus); 
    … …  
    329471 
    330472    // some double precision versions of the above values 
    331     df_A = sqrt(2) * 3.1415926535897931 * sqrt(n - 1.0/24.0); 
    332     df_B = 2.0 * sqrt(3) * (n - 1.0/24.0); 
    333     df_C = sqrt(2) * 3.1415926535897931 * sqrt(n - 1.0/24.0) / sqrt(3); 
    334     df_D = 2.0 * (n - 1.0/24.0) * sqrt(n - 1.0/24.0); 
     473 
     474    d_A = sqrt(2) * 3.1415926535897931 * sqrt(n - 1.0/24.0); 
     475    d_B = 2.0 * sqrt(3) * (n - 1.0/24.0); 
     476    d_C = sqrt(2) * 3.1415926535897931 * sqrt(n - 1.0/24.0) / sqrt(3); 
     477    d_D = 2.0 * (n - 1.0/24.0) * sqrt(n - 1.0/24.0); 
    335478 
    336479} 
    337480 
    338481void mp_f(mpfr_t result, unsigned int k) { 
     482    // compute f_n(k) as described in the introduction 
     483    // 
     484    // notice that this doesn't use n - the "constants" 
     485    // A, B, C, and D depend on n, but they are precomputed 
     486    // once before this function is called. 
     487     
    339488    //result =  pi * sqrt(2) * cosh(A/(sqrt(3)*k))/(B*k) - sinh(C/k)/D; 
    340     mpfr_set(result, mp_pi, round_mode);                    // result = pi 
    341     mpfr_mul(result, result, mp_sqrt2, round_mode);         // result = sqrt(2) * pi 
    342      
    343     mpfr_div(tempf1, mp_A, mp_sqrt3, round_mode);            // temp1 = mp_A/sqrt(3) 
    344     mpfr_div_ui(tempf1, tempf1, k, round_mode);               // temp1 = mp_A/(sqrt(3) * k) 
    345     mpfr_cosh(tempf1, tempf1, round_mode);                    // temp1 = cosh(mp_A/(sqrt(3) * k)) 
    346     mpfr_div(tempf1, tempf1, mp_B, round_mode);               // temp1 = cosh(mp_A/(sqrt(3) * k))/mp_B 
    347     mpfr_div_ui(tempf1, tempf1, k, round_mode);               // temp1 = cosh(mp_A/(sqrt(3) * k))/(mp_B*k) 
    348  
    349     mpfr_mul(result, result, tempf1, round_mode);            // result = sqrt(2) * pi * cosh(mp_A/(sqrt(3) * k))/(mp_B*k) 
    350      
    351     mpfr_div_ui(tempf1, mp_C, k, round_mode);                // temp1 = mp_C/k 
    352     mpfr_sinh(tempf1, tempf1, round_mode);                    // temp1 = sinh(mp_C/k) 
    353     mpfr_div(tempf1, tempf1, mp_D, round_mode);               // temp1 = sinh(mp_C/k)/D 
    354  
    355     mpfr_sub(result, result, tempf1, round_mode);            // result = RESULT! 
    356     //return 3.1415926535897931 * sqrt(2) * cosh(df_A/(sqrt(3)*k))/(df_B*k) - sinh(df_C/k)/df_D; 
    357 } 
    358  
    359 void mp_s(mpfr_t result, unsigned int h,unsigned int q) { 
    360     if(q < 3) { 
     489 
     490                                                                    // 
     491    mpfr_set(result, mp_pi, round_mode);                            // result = pi 
     492     
     493    mpfr_mul(result, result, mp_sqrt2, round_mode);                 // result = sqrt(2) * pi 
     494                                                                    // 
     495    mpfr_div(tempf1, mp_A, mp_sqrt3, round_mode);                   // temp1 = mp_A/sqrt(3) 
     496    mpfr_div_ui(tempf1, tempf1, k, round_mode);                     // temp1 = mp_A/(sqrt(3) * k) 
     497    mpfr_cosh(tempf1, tempf1, round_mode);                          // temp1 = cosh(mp_A/(sqrt(3) * k)) 
     498    mpfr_div(tempf1, tempf1, mp_B, round_mode);                     // temp1 = cosh(mp_A/(sqrt(3) * k))/mp_B 
     499    mpfr_div_ui(tempf1, tempf1, k, round_mode);                     // temp1 = cosh(mp_A/(sqrt(3) * k))/(mp_B*k) 
     500                                                                    // 
     501    mpfr_mul(result, result, tempf1, round_mode);                   // result = sqrt(2) * pi * cosh(mp_A/(sqrt(3) * k))/(mp_B*k) 
     502                                                                    // 
     503    mpfr_div_ui(tempf1, mp_C, k, round_mode);                       // temp1 = mp_C/k 
     504    mpfr_sinh(tempf1, tempf1, round_mode);                          // temp1 = sinh(mp_C/k) 
     505    mpfr_div(tempf1, tempf1, mp_D, round_mode);                     // temp1 = sinh(mp_C/k)/D 
     506                                                                    // 
     507    mpfr_sub(result, result, tempf1, round_mode);                   // result = RESULT! 
     508                                                                    // 
     509    //return pi * sqrt(2) * cosh(A/(sqrt(3)*k))/(B*k) - sinh(C/k)/D; 
     510} 
     511 
     512// call the following when 4k < sqrt(UINT_MAX) 
     513// 
     514// TODO: maybe a faster version of this can be written without using mpq_t, 
     515//       or maybe this version can be smarter. 
     516// 
     517//       The actual value of the cosine that we compute using s(h,k) 
     518//       only depends on {s(h,k)/2}, that is, the fractional 
     519//       part of s(h,k)/2. It may be possible to make use of this somehow. 
     520// 
     521//       NOTE: when we compute p(1000000000), 
     522//       it takes about 3m 30s (on my laptop). If we uncomment 
     523//       the first two lines below, so that this function doesn't actually 
     524//       compute anything, it takes about 3m 10s to run. So not that much time 
     525//       is spent in this function when computing for large n 
     526void q_s(mpq_t result, unsigned int h, unsigned int k) { 
     527    //mpq_set_ui(result, 1, 1); 
     528    //return; 
     529     
     530    if(k < 3) { 
     531        mpq_set_ui(result, 0, 1); 
     532        return; 
     533    } 
     534 
     535    if (h == 1) { 
     536        unsigned int d = GCD( (k-1)*(k-2), 12*k); 
     537        if(d > 1) { 
     538            mpq_set_ui(result, ((k-1)*(k-2))/d, (12*k)/d); 
     539        } 
     540        else { 
     541            mpq_set_ui(result, (k-1)*(k-2), 12*k); 
     542        } 
     543        //mpq_canonicalize(result); 
     544        return; 
     545    } 
     546    // TODO: In the function mp_s() there are a few special cases for special forms of h and k. 
     547    // (And there are more special cases listed in one of the references listed in the introduction.) 
     548    // 
     549    // It may be advantageous to implement some here, but I'm not sure 
     550    // if there is any real speed benefit to this. 
     551    // 
     552    // In the mpfr_t version of this function, the speedups didn't seem to help too much, but 
     553    // they might make more of a difference when using mpq_t. 
     554 
     555    // if h = 2 and k is odd, we have 
     556    // s(h,k) = (k-1)*(k-5)/24k 
     557    //if(h == 2 && k > 5 && k % 2 == 1) { 
     558    //    unsigned int d = GCD( (k-1)*(k-5), 24*k); 
     559    //    if(d > 1) { 
     560    //        mpq_set_ui(result, ((k-1)*(k-5))/d, (24*k)/d); 
     561    //    } 
     562    //    else { 
     563    //        mpq_set_ui(result, (k-1)*(k-5), 24*k); 
     564    //    } 
     565    //    return; 
     566    //} 
     567 
     568/*     
     569 
     570    // if k % h == 1, then 
     571    // 
     572    //      s(h,k) = (k-1)(k - h^2 - 1)/(12hk) 
     573    // 
     574 
     575    // We need to be a little careful here because k - h^2 - 1 can be negative. 
     576    if(k % h == 1) { 
     577        int num = (k-1)*(k - h*h - 1); 
     578        int den = 12*k*h; 
     579        int d = GCD(num, den); 
     580        if(d > 1) { 
     581            mpq_set_si(result, num/d, den/d); 
     582        } 
     583        else { 
     584            mpq_set_si(result, num, den); 
     585        } 
     586        return; 
     587    } 
     588 
     589    // if k % h == 2, then 
     590    // 
     591    //      s(h,k) = (k-2)[k - .5(h^2 + 1)]/(12hk) 
     592    // 
     593 
     594    //if(k % h == 2) { 
     595    //} 
     596*/ 
     597     
     598 
     599 
     600     
     601    mpq_set_ui(result, 0, 1);                             // result = 0 
     602 
     603    int r1 = k; 
     604    int r2 = h; 
     605 
     606    int n = 0; 
     607    int temp3; 
     608    while(r1 > 0 && r2 > 0) { 
     609        unsigned int d = GCD(r1 * r1 + r2 * r2 + 1, r1 * r2); 
     610        if(d > 1) { 
     611            mpq_set_ui(qtemp1, (r1 * r1 + r2 * r2 + 1)/d, (r1 * r2)/d); 
     612        } 
     613        else{ 
     614            mpq_set_ui(qtemp1, r1 * r1 + r2 * r2 + 1, r1 * r2); 
     615        } 
     616        //mpq_canonicalize(qtemp1); 
     617         
     618        if(n % 2 == 0){                                             // 
     619            mpq_add(result, result, qtemp1);                        // result += temp1; 
     620        }                                                           // 
     621        else {                                                      // 
     622            mpq_sub(result, result, qtemp1);                        // result -= temp1; 
     623        }                                                           // 
     624        temp3 = r1 % r2;                                            // 
     625        r1 = r2;                                                    // 
     626        r2 = temp3;                                                 // 
     627        n++;                                                        // 
     628    } 
     629 
     630    mpq_set_ui(qtemp1, 1, 12); 
     631    mpq_mul(result, result, qtemp1);                                // result = result * 1.0/12.0; 
     632     
     633     
     634    if(n % 2 == 1) { 
     635        mpq_set_ui(qtemp1, 1, 4); 
     636        mpq_sub(result, result, qtemp1);                            // result = result - .25; 
     637    } 
     638 
     639} 
     640 
     641 
     642void mp_s(mpfr_t result, unsigned int h,unsigned int k) { 
     643    // Compute s(h,k) using some clever tricks. 
     644 
     645    // If k < 3, then we know that the result is going to be 0. 
     646    if(k < 3) { 
    361647        mpfr_set_ui(result, 0, round_mode); 
    362648        return; 
    363649    } 
    364650 
    365     //double result, R1, R2, temp1, temp2; 
    366651    unsigned int n, r1, r2, temp3 = 0; 
     652    
     653    // If h = 1, then the result has the form 
     654    //  
     655    //      s(h,k) = (k-1)(k-2)/(12k). 
     656    // 
    367657     
    368658    if(h == 1) { 
    369         if(q < (unsigned int)(sqrt(UINT_MAX))) { 
    370             //elternate, assuming no overflow 
    371             mpfr_set_ui(result, (q-1)*(q-2), round_mode); 
    372             mpfr_div_ui(result, result, 12*q, round_mode); 
     659         
     660        // When k is small enough, we can do some of the computation using 
     661        // just ordinary integer arithmetic. 
     662         
     663        if(k < (unsigned int)(sqrt(UINT_MAX))) { 
     664            mpfr_set_ui(result, (k-1)*(k-2), round_mode); 
     665            mpfr_div_ui(result, result, 12*k, round_mode); 
    373666        } 
    374667        else { 
    375             mpz_set_ui(ztemp1, q-1);                                // temp = q-1 
    376             mpz_mul_ui(ztemp1, ztemp1, q-2);                        // temp = (q-1)(q-2) 
     668            // When k is very large, we might need to use this code. 
     669            mpz_set_ui(ztemp1, k-1);                                // temp = k-1 
     670            mpz_mul_ui(ztemp1, ztemp1, k-2);                        // temp = (k-1)(k-2) 
    377671         
    378             mpfr_set_z(result, ztemp1, round_mode);                 // result = (q-1)(q-2) 
    379             mpz_set_ui(ztemp1, q);                                  // temp = q 
    380             mpz_mul_ui(ztemp1, ztemp1, 12);                         // temp = 12q 
    381             mpfr_div_z(result, result, ztemp1, round_mode);         // result = (q-1)(q-2)/12q 
    382         } 
    383 /* 
    384         //result = (q-1)*(q-2)/(12*q); 
    385         mpfr_set_ui(result, q-1, round_mode); 
    386         mpfr_mul_ui(result, result, q-2, round_mode); 
    387         mpfr_div_ui(result, result, q, round_mode);         // in this step we don't want to assume that 12*q will not overflow 
    388         mpfr_div_ui(result, result, 12, round_mode); 
    389 */ 
    390  
    391  
     672            mpfr_set_z(result, ztemp1, round_mode);                 // result = (k-1)(k-2) 
     673            mpz_set_ui(ztemp1, k);                                  // temp = k 
     674            mpz_mul_ui(ztemp1, ztemp1, 12);                         // temp = 12k 
     675            mpfr_div_z(result, result, ztemp1, round_mode);         // result = (k-1)(k-2)/12k 
     676        } 
    392677        return; 
    393678    } 
    394679 
    395     mpfr_set_ui(result, 0, round_mode);             // result = 0 
    396  
    397     if(debugs) 
    398         if(h > q) { 
    399             cout << "oops in mp_s()" << endl; 
    400             exit(1); 
    401         } 
    402     r1 = q; 
     680    // TODO: Below are a few special cases for special forms of h and k. 
     681    // 
     682    // It may be advantageous to implement a few more, but I'm not even sure 
     683    // that there is any real speed benefit to the special forms that are here 
     684    // now. 
     685 
     686    // if h = 2 and k is odd, then s(h,k) is given by 
     687    // (we need k > 5 because k is unsigned. k = 3 or 5 should 
     688    // be special cased.) 
     689    // 
     690    //      s(h,k) = (k-1)(k-5)/(24k) 
     691    // 
     692    if(h == 2 && k > 5 && k % 2 == 1) { 
     693         
     694        // When k is small enough, we can do some of the computation using 
     695        // just ordinary integer arithmetic. 
     696         
     697        if(k < (unsigned int)(sqrt(UINT_MAX))) { 
     698            mpfr_set_ui(result, (k-1)*(k-5), round_mode); 
     699            mpfr_div_ui(result, result, 24*k, round_mode); 
     700        } 
     701        else { 
     702            // When k is very large, we might need to use this code. 
     703            mpz_set_ui(ztemp1, k-1);                                // temp = k-1 
     704            mpz_mul_ui(ztemp1, ztemp1, k-5);                        // temp = (k-1)(k-5) 
     705         
     706            mpfr_set_z(result, ztemp1, round_mode);                 // result = (k-1)(k-5) 
     707            mpz_set_ui(ztemp1, k);                                  // temp = k 
     708            mpz_mul_ui(ztemp1, ztemp1, 24);                         // temp = 24k 
     709            mpfr_div_z(result, result, ztemp1, round_mode);         // result = (k-1)(k-5)/24k 
     710        } 
     711        return; 
     712    } 
     713 
     714    // if k % h == 1, then 
     715    // 
     716    //      s(h,k) = (k-1)(k - h^2 - 1)/(12hk) 
     717    // 
     718 
     719    if(k % h == 1) { 
     720        if(4*k < (unsigned int)(sqrt(UINT_MAX))) { 
     721            mpfr_set_si(result, (k-1)*(k - h*h - 1), round_mode); 
     722            mpfr_div_ui(result, result, 12*h*k, round_mode); 
     723            return; 
     724        } 
     725    } 
     726 
     727    // if k % h == 2, then 
     728    // 
     729    //      s(h,k) = (k-2)[k - .5(h^2 + 1)]/(12hk) 
     730    // 
     731    // 
     732    // 
     733 
     734 
     735     
     736    // At this point we have given up hope of using a special form and fall back on our generic  
     737    // algorithm.  
     738 
     739 
     740    mpfr_set_ui(result, 0, round_mode);                             // result = 0 
     741 
     742    r1 = k; 
    403743    r2 = h; 
    404744 
    405745    n = 0; 
    406746    while(r1 > 0 && r2 > 0) { 
    407         if(r1 < (unsigned int)(sqrt(UINT_MAX)/2.0)) {       // if r1 is small enough we can use 
    408                                                         // standard C integers 
    409                                                         // NOTE: squareroot computation should be optimized by the compiler. 
     747        if(r1 < (unsigned int)(sqrt(UINT_MAX)/2.0)) {               // if r1 is small enough we can use 
     748                                                                    // standard C integers 
     749                                                                    // NOTE: squareroot computation should be optimized by the compiler. 
    410750            mpfr_set_ui(temps1, r1*r1 + r2*r2 + 1, round_mode); 
    411751            mpfr_div_ui(temps1, temps1, r1 * r2, round_mode); 
    … …  
    413753        } 
    414754        else { 
    415             //temp1 = (R1*R1 + R2*R2 + 1.0)/(R1 * R2);      // again, we are NOT going to assume no overflow 
    416  
    417             mpfr_set_ui(temps1, r1, round_mode);             // temp1 = r1 
    418             mpfr_mul_ui(temps1, temps1, r1, round_mode);      // temp1 = r1 * r1 
    419  
    420             mpfr_set_ui(temps2, r2, round_mode);             // temp2 = r2 
    421             mpfr_mul_ui(temps2, temps2, r2, round_mode);      // temp2 = r2 * r2 
    422  
    423             mpfr_add(temps1, temps1, temps2, round_mode);      // temp1 = r1*r1 + r2*r2 
    424             mpfr_add_ui(temps1, temps1, 1, round_mode);       // temp1 = r1*r1 + r2*r2 + 1 
    425  
    426             mpfr_div_ui(temps1, temps1, r1, round_mode);      // temp1 = (r1*r1 + r2*r2 + 1)/r1 
    427             mpfr_div_ui(temps1, temps1, r2, round_mode);      // temp1 = (r1*r1 + r2*r2 + 1)/(r1 * r2) 
    428         } 
    429         if(n % 2 == 0){ 
    430             mpfr_add(result, result, temps1, round_mode); // result += temp1; 
    431         } 
    432         else { 
    433             mpfr_sub(result, result, temps1, round_mode); // result -= temp1; 
    434         } 
    435         temp3 = r1 % r2; 
    436         r1 = r2; 
    437         r2 = temp3; 
    438         n++; 
    439     } 
    440  
    441     //cout << temps1 << endl; 
    442     //cout << result << endl; 
    443  
    444     mpfr_mul(result, result, mp_one_over_12, round_mode); // result = result * 1.0/12.0; 
     755            //temp1 = (R1*R1 + R2*R2 + 1.0)/(R1 * R2);              // 
     756                                                                    // 
     757            mpfr_set_ui(temps1, r1, round_mode);                    // temp1 = r1 
     758            mpfr_mul_ui(temps1, temps1, r1, round_mode);            // temp1 = r1 * r1 
     759                                                                    // 
     760            mpfr_set_ui(temps2, r2, round_mode);                    // temp2 = r2 
     761            mpfr_mul_ui(temps2, temps2, r2, round_mode);            // temp2 = r2 * r2 
     762                                                                    // 
     763            mpfr_add(temps1, temps1, temps2, round_mode);           // temp1 = r1*r1 + r2*r2 
     764            mpfr_add_ui(temps1, temps1, 1, round_mode);             // temp1 = r1*r1 + r2*r2 + 1 
     765                                                                    // 
     766            mpfr_div_ui(temps1, temps1, r1, round_mode);            // temp1 = (r1*r1 + r2*r2 + 1)/r1 
     767            mpfr_div_ui(temps1, temps1, r2, round_mode);            // temp1 = (r1*r1 + r2*r2 + 1)/(r1 * r2) 
     768        }                                                           // 
     769        if(n % 2 == 0){                                             // 
     770            mpfr_add(result, result, temps1, round_mode);           // result += temp1; 
     771        }                                                           // 
     772        else {                                                      // 
     773            mpfr_sub(result, result, temps1, round_mode);           // result -= temp1; 
     774        }                                                           // 
     775        temp3 = r1 % r2;                                            // 
     776        r1 = r2;                                                    // 
     777        r2 = temp3;                                                 // 
     778        n++;                                                        // 
     779    } 
     780 
     781 
     782    mpfr_mul(result, result, mp_one_over_12, round_mode);           // result = result * 1.0/12.0; 
    445783     
    446784     
    447785    if(n % 2 == 1) { 
    448786        mpfr_set_d(temps1, .25, round_mode); 
    449         mpfr_sub(result, result, temps1, round_mode);      // result = result - .25; 
    450     } 
    451  
    452  
    453  
    454     //return result; 
     787        mpfr_sub(result, result, temps1, round_mode);               // result = result - .25; 
     788    } 
     789 
    455790} 
    456791 
    457792 
    458793void mp_a(mpfr_t result, unsigned int n, unsigned int k) { 
     794    // compute a(n,k) 
    459795 
    460796    if (k == 1) { 
    461         mpfr_set_ui(result, 1, round_mode);     //result = 1 
     797        mpfr_set_ui(result, 1, round_mode);                         //result = 1 
    462798        return; 
    463799    } 
    … …  
    468804    for(h = 1; h < k+1; h++) { 
    469805        if(GCD(h,k) == 1) { 
    470             //result += cos(pi * ( s(h,k) - (2.0 * h * n)/k) ); 
    471             mp_s(tempa1, h, k);                              // temp1 = s(h,k) 
    472806             
    473             if(debugs) { 
    474                 cout << "s(" << h << "," << k << ") = "; 
    475                 mpfr_out_str(stdout, 10, 0, tempa1, round_mode); 
    476                 cout << endl; 
     807            // Note that we compute each term of the summand as 
     808            //      result += cos(pi * ( s(h,k) - (2.0 * h * n)/k) ); 
     809            // 
     810            // This is the real part of the exponential that was written 
     811            // down in the introduction, and we don't need to compute 
     812            // the imaginary part because we know that, in the end, the 
     813            // imaginary part will be 0, as we are computing an integer. 
     814             
     815            if(4*k < (unsigned int)(sqrt(UINT_MAX))) { 
     816                q_s(qtemp2, h, k); 
     817                 
     818                //mpfr_mul_q(tempa1, mp_pi, qtemp2, round_mode); 
     819                //mpfr_mul_ui(tempa1, tempa1, k * k, round_mode); 
     820 
     821                //mpfr_set_q(tempa1, qtemp2, round_mode); 
     822                unsigned int r = n % k;                                     // here we make use of the fact that the  
     823                unsigned int d = GCD(r,k);                                  // cos() term written above only depends 
     824                unsigned int K;                                             // on {hn/k}. 
     825                if(d > 1) { 
     826                    r = r/d; 
     827                    K = k/d; 
     828                } 
     829                else { 
     830                    K = k; 
     831                } 
     832                if(K % 2 == 0) { 
     833                    K = K/2; 
     834                } 
     835                else { 
     836                    r = r * 2; 
     837                } 
     838                mpq_set_ui(qtemp3, h*r, K); 
     839                mpq_sub(qtemp2, qtemp2, qtemp3); 
     840                /* 
     841                mpfr_set_q(tempa2, qtemp2, round_mode);                 // This might be faster, according to 
     842                cospi(tempa1, tempa2);                                  // the comments in Ralf Stephan's part.c, but 
     843                                                                        // I haven't noticed a significant speed up. 
     844                                                                        // (Perhaps a different version that takes an mpq_t 
     845                                                                        // as an input might be faster.) 
     846                */ 
     847                mpfr_mul_q(tempa1, mp_pi, qtemp2, round_mode); 
     848                mpfr_cos(tempa1, tempa1, round_mode); 
     849                mpfr_add(result, result, tempa1, round_mode); 
    477850            } 
    478  
    479             mpfr_set_ui(tempa2, 2, round_mode);              // temp2 = 2 
    480             mpfr_mul_ui(tempa2, tempa2, h, round_mode);       // temp2 = 2h 
    481             mpfr_mul_ui(tempa2, tempa2, n, round_mode);       // temp2 = 2hn 
    482             mpfr_div_ui(tempa2, tempa2, k, round_mode);       // temp2 = 2hn/k 
     851            else{ 
     852                mp_s(tempa1, h, k);                                     // temp1 = s(h,k) 
    483853             
    484             mpfr_sub(tempa1, tempa1, tempa2, round_mode);      // temp1 = s(h,k) - 2hn/k 
    485             mpfr_mul(tempa1, tempa1, mp_pi, round_mode);      // temp1 = pi * (s(h,k) - 2hn/k) 
    486             mpfr_cos(tempa1, tempa1, round_mode);             // temp1 = cos( ditto ) 
    487  
    488             mpfr_add(result, result, tempa1, round_mode);    // result = result + temp1 
    489         } 
     854                mpfr_set_ui(tempa2, 2, round_mode);                     // temp2 = 2 
     855                mpfr_mul_ui(tempa2, tempa2, h, round_mode);             // temp2 = 2h 
     856                mpfr_mul_ui(tempa2, tempa2, n, round_mode);             // temp2 = 2hn 
     857                mpfr_div_ui(tempa2, tempa2, k, round_mode);             // temp2 = 2hn/k 
     858             
     859                mpfr_sub(tempa1, tempa1, tempa2, round_mode);           // temp1 = s(h,k) - 2hn/k 
     860                mpfr_mul(tempa1, tempa1, mp_pi, round_mode);            // temp1 = pi * (s(h,k) - 2hn/k) 
     861                mpfr_cos(tempa1, tempa1, round_mode);                   // temp1 = cos( ditto ) 
     862 
     863                mpfr_add(result, result, tempa1, round_mode);           // result = result + temp1 
     864            } 
     865         
     866        } 
     867         
    490868    } 
    491869 
    … …  
    494872 
    495873void mp_t(mpfr_t result, unsigned int n) { 
     874    // This is the function that actually computes p(n). 
     875    // 
     876    // More specifically, it computes t(n,N) to within 1/2 of p(n), and then 
     877    // we can find p(n) by rounding. 
     878    // 
     879    // 
    496880    // NOTE: result should NOT have been initialized when this is called,  
    497881    // as we initialize it to the proper precision in this function. 
    498882 
    499     unsigned int initial_precision = compute_initial_precision(n); 
    500     mpfr_t t1, t2; 
    501     mpfr_init2(t1, initial_precision); 
    502     mpfr_init2(t2, initial_precision); 
    503     mpfr_init2(result, initial_precision); 
    504     mpfr_set_ui(result, 0, round_mode); 
    505  
    506     mp_init(initial_precision, n); 
    507     mp_set_precision(initial_precision); 
     883    unsigned int initial_precision = compute_initial_precision(n);  // We begin by computing the precision necessary to hold the final answer. 
     884                                                                    // and then initialize both the result and some temporary variables to that 
     885                                                                    // precision. 
     886    mpfr_t t1, t2;                                                  // 
     887    mpfr_init2(t1, initial_precision);                              // 
     888    mpfr_init2(t2, initial_precision);                              // 
     889    mpfr_init2(result, initial_precision);                          // 
     890    mpfr_set_ui(result, 0, round_mode);                             // 
     891                                                                     
     892    initialize_constants(initial_precision, n);                     // Now that we have the precision information, we initialize some constants 
     893                                                                    // that will be used throughout, and also 
     894                                                                    // 
     895    mp_set_precision(initial_precision);                            // set the precision of the "temp" variables that are used in individual functions. 
    508896 
    509897    unsigned int current_precision = initial_precision; 
    510898    unsigned int new_precision; 
    511899     
    512     double remainder = 100.0; 
    513 //    d_f_precompute(n); 
    514     unsigned int k = 1; 
    515     for(k = 1; current_precision > min_precision; k++) { 
     900    double remainder = 0.5772156649;                                // (We just need the remainder to be initialized to something. This 
     901                                                                    // seems like as good a number as any.) 
     902     
     903                                                                     
     904    // We start by computing with high precision arithmetic, until 
     905    // we are sure enough that we don't need that much precision 
     906    // anymore. Once current_precision == min_precision, we drop 
     907    // out of this loop and switch to a computation 
     908    // that only involves doubles. 
     909 
     910    unsigned int k = 1;                                             // (k holds the index of the summand that we are computing.) 
     911    for(k = 1; current_precision > min_precision; k++) {            // 
     912    /*    
     913        new_precision = compute_current_precision(n,k);             // After computing one summand, check what the new precision should be. 
     914        if(new_precision != current_precision) {                    // If the precision changes, we need to clear 
     915            current_precision = new_precision;                      // and reinitialize all "temp" variables to 
     916            mp_set_precision(current_precision);                    // use lower precision. 
     917            mpfr_clear(t1); mpfr_clear(t2);                         // 
     918            mpfr_init2(t1, current_precision);                      // 
     919            mpfr_init2(t2, current_precision);                      // 
     920        } 
     921     
     922        continue; 
     923    */ 
    516924        mpfr_sqrt_ui(t1, k, round_mode);                            // t1 = sqrt(k) 
    517          
     925                                                                    // 
    518926        mp_a(t2, n, k);                                             // t2 = A_k(n) 
    519927         
    520928        if(debuga) { 
    521929            cout << "a(" << k <<  ") = "; 
    522             mpfr_out_str(stdout, 10, 0, t2, round_mode); 
     930            mpfr_out_str(stdout, 10, 10, t2, round_mode); 
    523931            cout << endl; 
    524932        } 
    525933         
    526934        mpfr_mul(t1, t1, t2, round_mode);                           // t1 = sqrt(k)*A_k(n) 
    527  
     935                                                                    // 
    528936        mp_f(t2, k);                                                // t2 = f_k(n) 
    529937 
    … …  
    535943 
    536944        mpfr_mul(t1, t1, t2, round_mode);                           // t1 = sqrt(k)*A_k(n)*f_k(n) 
    537  
     945                                                                    // 
    538946        mpfr_add(result, result, t1, round_mode);                   // result += summand 
    539947 
    540948        if(debugt) { 
    541             cout << "Partial sum " << k << " = "; 
    542             mpfr_out_str(stdout, 10, 0, result, round_mode); 
     949            //cout << "Partial sum " << k << " = "; 
     950            //mpfr_out_str(stdout, 10, 0, result, round_mode); 
     951            //cout << endl; 
     952            int num_digits = 20; 
     953            int num_extra_digits; 
     954            char digits[num_digits + 1]; 
     955            num_extra_digits = grab_last_digits(digits, 5, t1); 
     956            grab_last_digits(digits, num_digits, result); 
     957 
     958            mpfr_out_str(stdout, 10, 10, t1, round_mode); 
    543959            cout << endl; 
    544         } 
    545         //K++; 
    546         //temp = sqrt(K) * d_A(n,k) * d_f(k); 
    547         //result += temp; 
    548  
    549         new_precision = compute_current_precision(n,k); 
    550         if(new_precision != current_precision) { 
    551             current_precision = new_precision; 
    552             mp_set_precision(current_precision); 
    553             mpfr_clear(t1); mpfr_clear(t2); 
    554             mpfr_init2(t1, current_precision); 
    555             mpfr_init2(t2, current_precision); 
    556         } 
    557     } 
    558  
    559     double result2 = 0; 
    560  
    561     for( ; remainder > .5; k++) { 
    562         result2 += sqrt(k) * d_A(n,k) * d_f(k); 
    563         remainder = compute_remainder(n,k); 
    564     } 
    565  
    566     mpfr_set_d(t1, result2, round_mode); 
    567     mpfr_add(result, result, t1, round_mode); 
    568  
    569 //    result = result/(3.1415926535897931 * sqrt(2)); 
    570     mpfr_div(result, result, mp_pi, round_mode); 
    571     mpfr_div(result, result, mp_sqrt2, round_mode); 
    572 } 
    573  
    574 //  double versions of the functions 
    575  
    576 void d_f_precompute(unsigned int n) { 
    577     df_A = sqrt(2) * 3.1415926535897931 * sqrt(n - 1.0/24.0); 
    578     df_B = 2.0 * sqrt(3) * (n - 1.0/24.0); 
    579     df_C = sqrt(2) * 3.1415926535897931 * sqrt(n - 1.0/24.0) / sqrt(3); 
    580     df_D = 2.0 * (n - 1.0/24.0) * sqrt(n - 1.0/24.0); 
    581 } 
     960 
     961            cout << k << ": current precision:"  << current_precision << ". 20 last digits of partial result: " << digits << ". Number of extra digits: " << num_extra_digits << endl; 
     962            cout.flush(); 
     963 
     964        } 
     965 
     966        new_precision = compute_current_precision(n,k);             // After computing one summand, check what the new precision should be. 
     967        if(new_precision != current_precision) {                    // If the precision changes, we need to clear 
     968            current_precision = new_precision;                      // and reinitialize all "temp" variables to 
     969            mp_set_precision(current_precision);                    // use lower precision. 
     970            mpfr_clear(t1); mpfr_clear(t2);                         // 
     971            mpfr_init2(t1, current_precision);                      // 
     972            mpfr_init2(t2, current_precision);                      // 
     973        } 
     974    } 
     975 
     976    double result2 = 0;                                             // (result2 will hold the result of the "tail end" 
     977                                                                    // computation using doubles.) 
     978 
     979    for( ; remainder > .5; k++) {                                   // (don't change k -- it is already the right value) 
     980        result2 += sqrt(k) * d_a(n,k) * d_f(k);                     // 
     981        remainder = compute_remainder(n,k);                         // Now we start computing the size of the remainder. Once 
     982                                                                    // it is small enough, we know that we have the answer. 
     983    } 
     984 
     985    mpfr_set_d(t1, result2, round_mode);                            // 
     986    mpfr_add(result, result, t1, round_mode);                       // We add together the main result and the tail end. 
     987 
     988    mpfr_div(result, result, mp_pi, round_mode);                    // The actual result is the sum that we have computed 
     989    mpfr_div(result, result, mp_sqrt2, round_mode);                 // divided by pi*sqrt(2). 
     990} 
     991 
     992//  Double versions of the functions, see the above functions for documentation. 
    582993 
    583994double d_f(unsigned int k) { 
    584     return 3.1415926535897931 * sqrt(2) * cosh(df_A/(sqrt(3)*k))/(df_B*k) - sinh(df_C/k)/df_D; 
    585 } 
    586  
    587 double d_s(unsigned int h,unsigned int q) { 
    588     if(q < 3) { 
     995    return 3.1415926535897931 * sqrt(2) * cosh(d_A/(sqrt(3)*k))/(d_B*k) - sinh(d_C/k)/d_D; 
     996} 
     997 
     998double d_s(unsigned int h,unsigned int k) { 
     999    if(k < 3) { 
    5891000        return 0.0; 
    5901001    } 
    … …  
    5941005     
    5951006    if(h == 1) { 
    596         double Q; 
    597         Q = q; 
    598         result = (Q-1)*(Q-2)/(12*Q); 
     1007        double K; 
     1008        K = k; 
     1009        result = (K-1)*(K-2)/(12*K); 
    5991010        return result; 
    6001011    } 
    6011012 
    6021013    result = 0; 
    603     R1 = q; 
     1014    R1 = k; 
    6041015    R2 = h; 
    6051016 
    606     r1 = q; 
     1017    r1 = k; 
    6071018    r2 = h; 
    6081019 
    … …  
    6331044} 
    6341045 
    635 double d_A(unsigned int n, unsigned int k) { 
     1046double d_a(unsigned int n, unsigned int k) { 
    6361047    double result; 
    6371048    result = 0; 
    … …  
    6831094 
    6841095 
     1096 
     1097 
     1098void cospi (mpfr_t res,  
     1099                mpfr_t x) 
     1100{ 
     1101//      mpfr_t t, tt, half, fourth; 
     1102         
     1103//      mpfr_init2 (t, prec); 
     1104//      mpfr_init2 (tt, prec); 
     1105//      mpfr_init2 (half, prec); 
     1106//      mpfr_init2 (fourth, prec); 
     1107 
     1108//      mpfr_set_ui (half, 1, r); 
     1109//      mpfr_div_2ui (half, half, 1, r); 
     1110//      mpfr_div_2ui (fourth, half, 1, r); 
     1111 
     1112        // NOTE: switched t to tempc2 
     1113        //       and tt to tempc1 
     1114 
     1115 
     1116        mp_rnd_t r = round_mode; 
     1117         
     1118 
     1119        mpfr_div_2ui (tempc1, x, 1, r); 
     1120        mpfr_floor (tempc1, tempc1); 
     1121        mpfr_mul_2ui (tempc1, tempc1, 1, r); 
     1122        mpfr_sub (tempc2, x, tempc1, r); 
     1123        if (mpfr_cmp_ui (tempc2, 1) > 0) 
     1124                mpfr_sub_ui (tempc2, tempc2, 2, r); 
     1125        mpfr_abs (tempc1, tempc2, r); 
     1126        if (mpfr_cmp (tempc1, half) > 0) 
     1127        { 
     1128                mpfr_ui_sub (tempc2, 1, tempc1, r); 
     1129                mpfr_abs (tempc1, tempc2, r); 
     1130                if (mpfr_cmp (tempc1, fourth) > 0) 
     1131                { 
     1132                        if (mpfr_sgn (tempc2) > 0) 
     1133                                mpfr_sub (tempc2, half, tempc2, r); 
     1134                        else 
     1135                                mpfr_add (tempc2, tempc2, half, r); 
     1136                        mpfr_mul (tempc2, tempc2, mp_pi, r); 
     1137                        mpfr_sin (tempc2, tempc2, r); 
     1138                } 
     1139                else 
     1140                { 
     1141                        mpfr_mul (tempc2, tempc2, mp_pi, r); 
     1142                        mpfr_cos (tempc2, tempc2, r);    
     1143                } 
     1144                mpfr_neg (res, tempc2, r); 
     1145        } 
     1146        else 
     1147        { 
     1148                mpfr_abs (tempc1, tempc2, r); 
     1149                if (mpfr_cmp (tempc1, fourth) > 0) 
     1150                { 
     1151                        if (mpfr_sgn (tempc2) > 0) 
     1152                                mpfr_sub (tempc2, half, tempc2, r); 
     1153                        else 
     1154                                mpfr_add (tempc2, tempc2, half, r); 
     1155                        mpfr_mul (tempc2, tempc2, mp_pi, r); 
     1156                        mpfr_sin (res, tempc2, r); 
     1157                } 
     1158                else 
     1159                { 
     1160                        mpfr_mul (tempc2, tempc2, mp_pi, r); 
     1161                        mpfr_cos (res, tempc2, r); 
     1162                } 
     1163        } 
     1164 
     1165//      mpfr_clear (half); 
     1166//      mpfr_clear (fourth); 
     1167//      mpfr_clear (t); 
     1168//      mpfr_clear (tt); 
     1169} 
     1170 
     1171 
     1172 
     1173void mpz_part(mpz_t result, unsigned int n) { 
     1174    if(n == 1) { 
     1175        mpz_set_str(result, "1", 10); 
     1176        return; 
     1177    } 
     1178     
     1179    mpfr_t mp_result; 
     1180 
     1181    mp_t(mp_result, n); 
     1182     
     1183    mpfr_get_z(result, mp_result, round_mode); 
     1184 
     1185    mpfr_clear(mp_result); 
     1186     
     1187    return; 
     1188} 
     1189 
     1190 
     1191 
     1192int main(int argc, char *argv[]){ 
     1193    //init(); 
     1194     
     1195    unsigned int n = 10; 
     1196 
     1197    if(argc > 1) 
     1198        n = atoi(argv[1]); 
     1199    else { 
     1200        cout << test() << endl; 
     1201        return 0; 
     1202    } 
     1203    //mpfr_t result; 
     1204     
     1205    //mp_t(result, n); 
     1206     
     1207    mpz_t answer; 
     1208    mpz_init(answer); 
     1209    mpz_part(answer, n); 
     1210 
     1211    //mpfr_get_z(answer, result, round_mode); 
     1212 
     1213    mpz_out_str (stdout, 10, answer); 
     1214     
     1215    cout << endl; 
     1216 
     1217    return 0; 
     1218} 
     1219 
     1220 
     1221bool test(bool longtest) { 
     1222    // The values given below are confirmed by multiple sources, so are probably correct. 
     1223    // TODO: There should be some more code here to test that answers satisfy the proper congruences that the should 
     1224    //  satisfy. On the other hand, it might be better to test this file from within SAGE, 
     1225    //  since that might be easier, and would certainly be more in line with the rest of 
     1226    //  SAGE. 
     1227 
     1228    mpz_t expected_value; 
     1229    mpz_t actual_value; 
     1230     
     1231    mpz_init(expected_value); 
     1232    mpz_init(actual_value); 
     1233 
     1234    // n = 1 
     1235    mpz_set_str(expected_value, "1", 10); 
     1236    mpz_part(actual_value, 1); 
     1237 
     1238    if(mpz_cmp(expected_value, actual_value) != 0) 
     1239        return false; 
     1240 
     1241    // n = 10 
     1242    mpz_set_str(expected_value, "42", 10); 
     1243    mpz_part(actual_value, 10); 
     1244 
     1245    if(mpz_cmp(expected_value, actual_value) != 0) 
     1246        return false; 
     1247     
     1248    // n = 100 
     1249    mpz_set_str(expected_value, "190569292", 10); 
     1250    mpz_part(actual_value, 100); 
     1251 
     1252    if(mpz_cmp(expected_value, actual_value) != 0) 
     1253        return false; 
     1254 
     1255    // n = 1000 
     1256    mpz_set_str(expected_value, "24061467864032622473692149727991", 10); 
     1257    mpz_part(actual_value, 1000); 
     1258 
     1259    if(mpz_cmp(expected_value, actual_value) != 0) 
     1260        return false; 
     1261 
     1262    // n = 10000 
     1263    mpz_set_str(expected_value, "36167251325636293988820471890953695495016030339315650422081868605887952568754066420592310556052906916435144", 10); 
     1264    mpz_part(actual_value, 10000); 
     1265 
     1266    if(mpz_cmp(expected_value, actual_value) != 0) 
     1267        return false; 
     1268 
     1269    // n = 100000 
     1270    mpz_set_str(expected_value, "27493510569775696512677516320986352688173429315980054758203125984302147328114964173055050741660736621590157844774296248940493063070200461792764493033510116079342457190155718943509725312466108452006369558934464248716828789832182345009262853831404597021307130674510624419227311238999702284408609370935531629697851569569892196108480158600569421098519", 10); 
     1271    mpz_part(actual_value, 100000); 
     1272 
     1273    if(mpz_cmp(expected_value, actual_value) != 0) 
     1274        return false; 
     1275     
     1276    // n = 1000000 
     1277    mpz_set_str(expected_value, "1471684986358223398631004760609895943484030484439142125334612747351666117418918618276330148873983597555842015374130600288095929387347128232270327849578001932784396072064228659048713020170971840761025676479860846908142829356706929785991290519899445490672219997823452874982974022288229850136767566294781887494687879003824699988197729200632068668735996662273816798266213482417208446631027428001918132198177180646511234542595026728424452592296781193448139994664730105742564359154794989181485285351370551399476719981691459022015599101959601417474075715430750022184895815209339012481734469448319323280150665384042994054179587751761294916248142479998802936507195257074485047571662771763903391442495113823298195263008336489826045837712202455304996382144601028531832004519046591968302787537418118486000612016852593542741980215046267245473237321845833427512524227465399130174076941280847400831542217999286071108336303316298289102444649696805395416791875480010852636774022023128467646919775022348562520747741843343657801534130704761975530375169707999287040285677841619347472368171772154046664303121315630003467104673818", 10); 
     1278    mpz_part(actual_value, 1000000); 
     1279 
     1280    if(mpz_cmp(expected_value, actual_value) != 0) 
     1281        return false; 
     1282     
     1283 
     1284 
     1285 
     1286    mpz_clear(expected_value); 
     1287    mpz_clear(actual_value); 
     1288 
     1289    return true; 
     1290} 
     1291 
     1292 
    6851293/* answer must have already been mpz_init'd. */ 
    6861294int part(mpz_t answer, unsigned int n){ 
Note: See TracChangeset for help on using the changeset viewer.