/* * PI-1.C * * DEMO 'C' PROGRAM TO COMPUTE PI USING * BORWEIN'S QUARTICALLY CONVERGENT ALGORITHM * * * Output from 'pi-1 200' : (I line wrapped for this page). * * start Borwein-Quartic PI calculation * PI now known to 2 digits of accuracy ... * PI now known to 8 digits of accuracy ... * PI now known to 18 digits of accuracy ... * PI now known to 38 digits of accuracy ... * PI now known to 82 digits of accuracy ... * PI now known to 168 digits of accuracy ... * Borwein-Quartic PI calculation complete * * Borwein PI = [3.1415926535897932384626433832795028841971693993751 * 058209749445923078164062862089986280348253421170679 * 821480865132823066470938446095505822317253594081284 * 8111745028410270193852110555964462294895493038196E+0] */ #include#include #include #include "m_apm.h" #define FALSE 0 #define TRUE 1 extern void compute_PI(M_APM, int); int main(argc, argv) int argc; char *argv[]; { int ii, dplaces; char *buffer; M_APM pi_mapm; /* declare the M_APM variable ... */ if (argc < 2) { fprintf(stdout,"\nUsage: pi-1 digits\t\t\t\t[Version 1.1]\n"); fprintf(stdout, " Compute PI to \'digits\' decimal places of accuracy\n"); exit(4); } pi_mapm = m_apm_init(); /* now initialize the M_APM variable ... */ dplaces = atoi(argv[1]); if (dplaces < 2) { fprintf(stdout,"\nUsage: pi-1 digits \n"); exit(6); } fprintf(stderr, "\nstart Borwein-Quartic PI calculation \n"); compute_PI(pi_mapm, dplaces); fprintf(stderr, "Borwein-Quartic PI calculation complete \n\n"); /* * find out how many digits so we can * convert the entire value into a string. * (plus some pad for decimal point, exponent, etc) */ ii = m_apm_significant_digits(pi_mapm) + 12; if ((buffer = (char *)malloc(ii)) == NULL) { fprintf(stdout,"PI demo: Out of memory \n"); exit(8); } /* now convert the MAPM number into a string */ m_apm_to_string(buffer, dplaces, pi_mapm); printf("Borwein PI = [%s] \n",buffer); m_apm_free(pi_mapm); free(buffer); exit(0); } /****************************************************************************/ /* * Calculate PI using the Borwein's Quartically Convergent Algorithm * * Init : U0 = sqrt(2) * B0 = 0 * P0 = 2 + sqrt(2) * * Iterate: * * * U = 0.5 * [ sqrt(U ) + 1 / sqrt(U ) ] * k+1 k k * * * sqrt(U ) * [ 1 + B ] * k k * B = ----------------------- * k+1 U + B * k k * * * P * B * [ 1 + U ] * k k+1 k+1 * P = ---------------------------- * k+1 1 + B * k+1 * * * P ---> PI * k+1 * */ void compute_PI(outv, places) M_APM outv; int places; { M_APM tmp1, tmp2, tmp3, c0_5, u0, u1, b0, b1, p0, p1; int ct, nn, dplaces, dflag; tmp1 = m_apm_init(); tmp2 = m_apm_init(); tmp3 = m_apm_init(); c0_5 = m_apm_init(); u0 = m_apm_init(); u1 = m_apm_init(); b0 = m_apm_init(); b1 = m_apm_init(); p0 = m_apm_init(); p1 = m_apm_init(); dplaces = places + 16; /* compute a few extra digits */ /* in the intermediate math */ dflag = FALSE; ct = 0; m_apm_set_string(c0_5, "0.5"); /* need this constant */ m_apm_sqrt(u0, dplaces, MM_Two); m_apm_add(p0, u0, MM_Two); /* b0 initialized to 0 from m_apm_init() */ while (TRUE) { m_apm_sqrt(tmp1, dplaces, u0); m_apm_reciprocal(tmp2, dplaces, tmp1); m_apm_add(tmp3, tmp1, tmp2); m_apm_multiply(u1, c0_5, tmp3); m_apm_add(tmp2, MM_One, b0); m_apm_multiply(tmp3, tmp1, tmp2); m_apm_add(tmp2, b0, u0); m_apm_divide(b1, dplaces, tmp3, tmp2); m_apm_add(tmp3, MM_One, u1); m_apm_multiply(tmp2, b1, tmp3); m_apm_multiply(tmp1, p0, tmp2); m_apm_add(tmp3, MM_One, b1); m_apm_divide(p1, dplaces, tmp1, tmp3); if (dflag) break; /* * compute the difference from this * iteration to the last one. */ m_apm_subtract(tmp2, p0, p1); if (++ct >= 4) { /* * if diff == 0, we're done. */ if (m_apm_sign(tmp2) == 0) break; } /* * if the exponent of the error term (small error like 2.47...E-65) * is small enough, break out after the *next* p1 is calculated. * * get the exponent of the error term, which will be negative. */ nn = -m_apm_exponent(tmp2); /* * normally, this wouldn't be here. it's nice in the demo though. */ fprintf(stderr, "PI now known to %d digits of accuracy ... \n",(2 * nn)); if ((4 * nn) >= dplaces) dflag = TRUE; /* set up for the next iteration */ m_apm_copy(b0, b1); m_apm_copy(u0, u1); m_apm_copy(p0, p1); } /* round to the accuracy asked for */ m_apm_round(outv, places, p1); m_apm_free(p1); m_apm_free(p0); m_apm_free(b1); m_apm_free(b0); m_apm_free(u1); m_apm_free(u0); m_apm_free(c0_5); m_apm_free(tmp3); m_apm_free(tmp2); m_apm_free(tmp1); } /****************************************************************************/ /****************************************************************************/ /* * PI-2.CPP * * DEMO 'C++' PROGRAM TO COMPUTE PI USING * BORWEIN'S QUARTICALLY CONVERGENT ALGORITHM * * * Output from 'pi-2 200' : (I line wrapped for this page). * * start Borwein-Quartic PI calculation * PI now known to 2 digits of accuracy ... * PI now known to 8 digits of accuracy ... * PI now known to 18 digits of accuracy ... * PI now known to 38 digits of accuracy ... * PI now known to 82 digits of accuracy ... * PI now known to 168 digits of accuracy ... * Borwein-Quartic PI calculation complete * * Borwein PI = [3.1415926535897932384626433832795028841971693993751 * 058209749445923078164062862089986280348253421170679 * 821480865132823066470938446095505822317253594081284 * 8111745028410270193852110555964462294895493038196E+0] */ #include #include #include #include "m_apm.h" #define FALSE 0 #define TRUE 1 extern MAPM compute_PI(int); int main(int argc, char **argv) { int ii, dplaces; char *buffer; MAPM pi_mapm; /* declare the MAPM variable ... */ if (argc < 2) { fprintf(stdout,"\nUsage: pi-2 digits\t\t\t\t[Version 1.1]\n"); fprintf(stdout, " Compute PI to \'digits\' decimal places of accuracy\n"); exit(4); } dplaces = atoi(argv[1]); if (dplaces < 2) { fprintf(stdout,"\nUsage: pi-2 digits \n"); exit(6); } fprintf(stderr, "\nstart Borwein-Quartic PI calculation \n"); pi_mapm = compute_PI(dplaces); fprintf(stderr, "Borwein-Quartic PI calculation complete \n\n"); /* * find out how many digits so we can * convert the entire value into a string. * (plus some pad for decimal point, exponent, etc) */ ii = pi_mapm.significant_digits() + 12; if ((buffer = (char *)malloc(ii)) == NULL) { fprintf(stdout,"PI demo: Out of memory \n"); exit(8); } /* now convert the MAPM number into a string */ pi_mapm.toString(buffer, dplaces); printf("Borwein PI = [%s] \n",buffer); free(buffer); exit(0); } /****************************************************************************/ /* * Calculate PI using the Borwein's Quartically Convergent Algorithm * * Init : U0 = sqrt(2) * B0 = 0 * P0 = 2 + sqrt(2) * * Iterate: * * * U = 0.5 * [ sqrt(U ) + 1 / sqrt(U ) ] * k+1 k k * * * sqrt(U ) * [ 1 + B ] * k k * B = ----------------------- * k+1 U + B * k k * * * P * B * [ 1 + U ] * k k+1 k+1 * P = ---------------------------- * k+1 1 + B * k+1 * * * P ---> PI * k+1 * */ MAPM compute_PI(int places) { MAPM t1, c0_5, u0, u1, b0, b1, p0, p1; int ct, nn, dplaces, dflag; dplaces = places + 16; // compute a few extra digits // in the intermediate math m_apm_cpp_precision(dplaces); // tell C++ about our precision dflag = FALSE; ct = 0; c0_5 = "0.5"; u0 = sqrt("2"); // make sure we use the MAPM sqrt b0 = 0; p0 = 2 + u0; while (TRUE) { /* * U = 0.5 * [ sqrt(U ) + 1 / sqrt(U ) ] * k+1 k k */ t1 = sqrt(u0); u1 = c0_5 * (t1 + 1 / t1); /* * sqrt(U ) * [ 1 + B ] * k k * B = ----------------------- * k+1 U + B * k k * */ b1 = (t1 * (1 + b0)) / (u0 + b0); /* * P * B * [ 1 + U ] * k k+1 k+1 * P = ---------------------------- * k+1 1 + B * k+1 */ p1 = (p0 * b1 * (1 + u1)) / (1 + b1); if (dflag) break; /* * compute the difference from this * iteration to the last one. */ t1 = p1 - p0; if (++ct >= 4) { /* * if diff == 0, we're done. */ if (t1 == 0) break; } /* * if the exponent of the error term (small error like 2.47...E-65) * is small enough, break out after the *next* p1 is calculated. * * get the exponent of the error term, which will be negative. */ nn = -(t1.exponent()); /* * normally, this wouldn't be here. it's nice in the demo though. */ fprintf(stderr, "PI now known to %d digits of accuracy ... \n",(2 * nn)); if ((4 * nn) >= dplaces) dflag = TRUE; /* set up for the next iteration */ /* */ /* also, keep the number of significant */ /* digits from growing without bound */ b0 = b1.round(dplaces); u0 = u1.round(dplaces); p0 = p1.round(dplaces); } /* round to the accuracy asked for */ return (p1.round(places)); } /****************************************************************************/ /****************************************************************************/ /* * PI-3.C * * DEMO 'C' PROGRAM TO COMPUTE PI USING * J.BORWEIN-P.BORWEIN'S ALGORITHM (1995) * */ #include #include #include #include "m_apm.h" #define FALSE 0 #define TRUE 1 #define Mult m_apm_multiply #define Add m_apm_add #define Subtract m_apm_subtract #define Divide m_apm_divide #define AssignString m_apm_set_string extern void compute_PI_2(M_APM, int); void m_apm_square(M_APM x, M_APM y) /* x = y^2 */ { m_apm_multiply(x, y, y); } /* * NOTE : The following 2 functions modify the * input value ('y'). Use with caution * in another application! */ /* x = y^4 */ void m_apm_pow4(M_APM x, M_APM y) { m_apm_square(x, y); m_apm_copy(y, x); m_apm_square(x, y); } /* x = sqrt(sqrt(y)) */ void m_apm_root4(M_APM x, int dec_places, M_APM y) { m_apm_sqrt(x, dec_places, y); m_apm_copy(y, x); m_apm_sqrt(x, dec_places, y); } int main(argc, argv) int argc; char *argv[]; { int ii, dplaces; char *buffer; M_APM pi_mapm; /* declare the M_APM variable ... */ if (argc < 2) { fprintf(stdout,"\nUsage: pi-3 digits\t\t\t\t[Version 1.1]\n"); fprintf(stdout, " Compute PI to \'digits\' decimal places of accuracy\n"); exit(4); } pi_mapm = m_apm_init(); /* now initialize the M_APM variable ... */ dplaces = atoi(argv[1]); if (dplaces < 2) { fprintf(stdout,"\nUsage: pi digits \n"); exit(6); } fprintf(stderr, "\nstart Borwein-Borwein PI calculation \n"); compute_PI_2(pi_mapm, dplaces); fprintf(stderr, "Borwein-Borwein PI calculation complete \n\n"); /* * find out how many digits so we can * convert the entire value into a string. * (plus some pad for decimal point, exponent, etc) */ ii = m_apm_significant_digits(pi_mapm) + 12; if ((buffer = (char *)malloc(ii)) == NULL) { fprintf(stdout,"PI demo: Out of memory \n"); exit(8); } /* now convert the MAPM number into a string */ m_apm_to_string(buffer, dplaces, pi_mapm); printf("Borwein PI = [%s] \n",buffer); m_apm_free(pi_mapm); free(buffer); exit(0); } /****************************************************************************/ /* * Calculate PI using the J. Borwein - P. Borwein's Algorithm (1995) * * * Init : A0 = 6 - 4 * sqrt(2) * B0 = sqrt(2) - 1 * * Iterate: * * 4 * 1 - quart( 1 - B ) * k * B = -------------------- * k+1 4 * 1 + quart( 1 - B ) * k * * * 4 2k+3 2 * A = A ( 1 + B ) - 2 B ( 1 + B + B ) * k+1 k k+1 k+1 k+1 k+1 * * * * * where quart(x) = sqrt(sqrt(x)) = 4th root of x * * 1 * A ---> ----- * k+1 PI * */ void compute_PI_2(outv, places) M_APM outv; int places; { M_APM tmp1, tmp2, tmp3, tmp4, a0, a1, b0, b1, MM_Six; int kk, ct, nn, dplaces, dflag; tmp1 = m_apm_init(); tmp2 = m_apm_init(); tmp3 = m_apm_init(); tmp4 = m_apm_init(); a0 = m_apm_init(); a1 = m_apm_init(); b0 = m_apm_init(); b1 = m_apm_init(); MM_Six = m_apm_init(); dplaces = places + 16; /* compute a few extra digits */ /* in the intermediate math */ dflag = FALSE; ct = 0; kk = 0; AssignString(MM_Six, "6"); m_apm_sqrt(tmp1, dplaces, MM_Two); Mult(tmp2, tmp1, MM_Four); Subtract(a0, MM_Six, tmp2); Subtract(b0, tmp1, MM_One); while (TRUE) { m_apm_pow4(tmp1, b0); /* tmp1 = b0^4 */ Subtract(tmp2, MM_One, tmp1); /* tmp2 = 1 - b0^4 */ m_apm_root4(tmp3, dplaces, tmp2); /* tmp3 = sqrt[4](1 - b0^4) */ Subtract(tmp1, MM_One, tmp3); /* tmp1 = 1 - tmp3 */ Add(tmp2, MM_One, tmp3); /* tmp2 = 1 + tmp3 */ Divide(b1, dplaces, tmp1, tmp2); /* b1 = tmp1 / tmp2 */ Add(tmp1, MM_One, b1); /* tmp1 = 1 + b1 */ m_apm_pow4(tmp2, tmp1); /* tmp2 = (1 + b1)^4 */ Mult(tmp3, a0, tmp2); /* tmp3 = a0 * (1 + b1)^4 */ m_apm_square(tmp1, b1); /* tmp1 = b1^2 */ Add(tmp2, tmp1, b1); /* tmp2 = b1 + b1^2 */ Add(tmp1, tmp2, MM_One); /* tmp1 = 1 + b1 + b1^2 */ Mult(tmp2, tmp1, b1); /* tmp2 = b1(1 + b1 + b1^2) */ m_apm_integer_pow(tmp1, dplaces, MM_Two, (2 * kk + 3)); /* tmp1 = 2^(2k+3) */ Mult(tmp4, tmp1, tmp2); /* tmp4 = 2^(2k+3) b1 (1 + b1 + b1^2) */ Subtract(tmp1, tmp3, tmp4); /* a1 = a0(1 + b1)^4 - tmp4 */ m_apm_round(a1, dplaces, tmp1); if (dflag) break; /* * compute the difference from this * iteration to the last one. */ m_apm_subtract(tmp2, a0, a1); if (++ct >= 4) { /* * if diff == 0, we're done. */ if (m_apm_sign(tmp2) == 0) break; } /* * if the exponent of the error term (small error like 2.47...E-65) * is small enough, break out after the *next* a1 is calculated. * * get the exponent of the error term, which will be negative. */ nn = -m_apm_exponent(tmp2); /* * normally, this wouldn't be here. it's nice in the demo though. */ fprintf(stderr, "PI now known to %d digits of accuracy ... \n",(4 * nn)); if ((15 * nn) >= dplaces) dflag = TRUE; /* set up for the next iteration */ m_apm_copy(b0, b1); m_apm_copy(a0, a1); kk++; } /* round to the accuracy asked for after taking the reciprocal */ m_apm_reciprocal(tmp1, dplaces, a1); m_apm_round(outv, places, tmp1); m_apm_free(MM_Six); m_apm_free(b1); m_apm_free(b0); m_apm_free(a1); m_apm_free(a0); m_apm_free(tmp4); m_apm_free(tmp3); m_apm_free(tmp2); m_apm_free(tmp1); } /****************************************************************************/ /****************************************************************************/ /* * PI-4.CPP * * DEMO 'C++' PROGRAM TO COMPUTE PI USING * J.BORWEIN-P.BORWEIN'S ALGORITHM (1995) * */ #include #include #include #include "m_apm.h" #define FALSE 0 #define TRUE 1 extern MAPM compute_PI_2(int); int main(int argc, char *argv[]) { int ii, dplaces; char *buffer; MAPM pi_mapm; /* declare the MAPM variable ... */ if (argc < 2) { fprintf(stdout,"\nUsage: pi-4 digits\t\t\t\t[Version 1.1]\n"); fprintf(stdout, " Compute PI to \'digits\' decimal places of accuracy\n"); exit(4); } dplaces = atoi(argv[1]); if (dplaces < 2) { fprintf(stdout,"\nUsage: pi-4 digits \n"); exit(6); } fprintf(stderr, "\nstart Borwein-Borwein PI calculation \n"); pi_mapm = compute_PI_2(dplaces); fprintf(stderr, "Borwein-Borwein PI calculation complete \n\n"); /* * find out how many digits so we can * convert the entire value into a string. * (plus some pad for decimal point, exponent, etc) */ ii = pi_mapm.significant_digits() + 12; if ((buffer = (char *)malloc(ii)) == NULL) { fprintf(stdout,"PI demo: Out of memory \n"); exit(8); } /* now convert the MAPM number into a string */ pi_mapm.toString(buffer, dplaces); printf("Borwein PI = [%s] \n",buffer); free(buffer); exit(0); } /****************************************************************************/ /* * Calculate PI using the J. Borwein - P. Borwein's Algorithm (1995) * * * Init : A0 = 6 - 4 * sqrt(2) * B0 = sqrt(2) - 1 * * Iterate: * * 4 * 1 - quart( 1 - B ) * k * B = -------------------- * k+1 4 * 1 + quart( 1 - B ) * k * * * 4 2k+3 2 * A = A ( 1 + B ) - 2 B ( 1 + B + B ) * k+1 k k+1 k+1 k+1 k+1 * * * * * where quart(x) = sqrt(sqrt(x)) = 4th root of x * * 1 * A ---> ----- * k+1 PI * */ MAPM compute_PI_2(int places) { MAPM t1, t2, c2, a0, a1, b0, b1; int kk, nn, dplaces, dflag; dplaces = places + 16; // compute a few extra digits // in the intermediate math m_apm_cpp_precision(dplaces); // tell C++ our precision requirements dflag = FALSE; kk = 0; t1 = sqrt("2.0"); // make sure we use the MAPM sqrt a0 = 6 - 4 * t1; b0 = t1 - 1; c2 = 2; while (TRUE) { t1 = b0.ipow(4); // t1 = b0 ^ 4 t1 = sqrt(sqrt(1 - t1)); // t1 = root_4(1 - t1) t1 = t1.round(dplaces); b1 = (1 - t1) / (1 + t1); b1 = b1.round(dplaces); t2 = 1 + b1; t2 = t2.ipow(4); // t2 = (1 + b1) ^ 4 t2 = t2.round(dplaces); t1 = c2.ipow(2 * kk + 3); // t1 = 2 ^ (2*kk+3) a1 = a0 * t2 - t1 * b1 * (1 + b1 + b1 * b1); a1 = a1.round(dplaces); if (dflag) break; // compute the difference from this // iteration to the last one. t1 = a1 - a0; if (kk >= 3) { // if diff == 0, we're done. if (t1 == 0) break; } /* * if the exponent of the error term (small error like 2.47...E-65) * is small enough, break out after the *next* a1 is calculated. * * get the exponent of the error term, which will be negative. */ nn = -(t1.exponent()); // normally, this wouldn't be here. it's nice in the demo though. fprintf(stderr, "PI now known to %d digits of accuracy ... \n",(4 * nn)); if ((15 * nn) >= dplaces) dflag = TRUE; // set up for the next iteration b0 = b1; a0 = a1; kk++; } // round to the accuracy asked for after taking the reciprocal t1 = 1 / a1; return(t1.round(places)); } /****************************************************************************/