Picarte's iteration
Experimental data
Claudio Gutiérrez, Mauricio Monsalve
Results
Sourcecode
<![CDATA[
/* iterations.c */
#include "stdio.h"
#include <sys/time.h>
#include <gmp.h>
/*Receives b and num, the number of iterations*/
/*It stores the result in inv. inv equals to the last x_i.*/
/*It returns the precision k.*/
int
picarte(mpz_t inv,mpz_t b,int num) {
mpz_t cajon,aux1,aux2,r,r2r;
mpz_init(cajon);
mpz_init(aux1);
mpz_init(aux2);
mpz_init(r);
mpz_init(r2r);
int i=1;
mpz_set_ui(inv,0);
mpz_set_ui(r,2);
//Iteracion de Picarte
for(;num>0;num--) {
mpz_mul_2exp(aux1,inv,i); /* 2^i x_i */
mpz_mul(aux2,inv,r); /* r_i x_i */
mpz_add(inv,aux1,aux2);
mpz_mul(r2r,r,r);
mpz_fdiv_qr(aux1,r,r2r,b); /* floor ( r_i^2/b ) */
mpz_add(aux2,inv,aux1);
mpz_swap(inv,aux2);
i=2*i;
//gmp_printf("%Zd %Zd %Zd \n",inv,r,aux1);
}
mpz_clear(cajon);
mpz_clear(aux1);
mpz_clear(aux2);
mpz_clear(r);
mpz_clear(r2r);
return i;
}
/*The same as Picarte*/
int
newton(mpz_t inv,mpz_t b,int num) {
mpz_t cajon,aux1,aux2,r,r2r,x2x;
mpz_init(cajon);
mpz_init(aux1);
mpz_init(aux2);
mpz_init(r);
mpz_init(r2r);
mpz_init(x2x);
int i=1;
mpz_set_ui(inv,0);
mpz_set_ui(r,2);
//Iteracion de Newton
for(;num>0;num--) {
mpz_mul_2exp(aux1,inv,i+1); /* 2^i 2 x_i */
mpz_mul(x2x,inv,inv); /* x_i^2 */
mpz_mul(aux2,x2x,b); /* b x_i^2 */
mpz_sub(inv,aux1,aux2);
mpz_mul(r2r,r,r);
mpz_fdiv_qr(aux1,r,r2r,b); /* floor ( r_i^2/b ) */
mpz_add(aux2,inv,aux1);
mpz_set(inv,aux2);
i=2*i;
//gmp_printf("%Zd %Zd %Zd \n",inv,r,aux1);
}
mpz_clear(cajon);
mpz_clear(aux1);
mpz_clear(aux2);
mpz_clear(r);
mpz_clear(r2r);
mpz_clear(x2x);
return i;
}
]]>