Prime numbers - new criteria

If we want to determine all the prime numbers smaller than n we have to check the divisors till sqrt{N-p}; here p is the first divisor of the greather number smaller than n.
In this case, for n=10000000(10mil), we make with 7901 steps lesser than the classic method - to go till sqrt{N}.

The classic criteria:

//21086848011 - the total number of steps
#include
#include
#include

unsigned long int i, n, p;
unsigned long int V;
long double S=0;

FILE * f;

void main(void)
{
//printf("Intr n:");
//scanf("%d",&n);

f=fopen("orig.txt","w+");

for(n=2;n<10000000;n++)
{
p=1;
V=(int)ceil(pow(n,0.5));
//printf("\nV=%d\n",V);
for(i=2;i<=V+1;i++)
{
if((n%i)==0) p=0;
S++;
}
if(p==0)
;//printf("NOK");
else
fprintf(f,"%Ld\n",n);
}
printf("\nAm terminat...");
fprintf(f,"S=%.1Lf",S);
fclose(f);
getch();

}

The new criteria:

//21086840110 - the total number of steps
//Diff==7901
#include
#include
#include

unsigned long int i, n, p;
unsigned long int V, X, G;
long double S=0;

FILE * f;

void main(void)
{
//printf("Intr n:");
//scanf("%d",&n);

f=fopen("prime.txt","w+");

for(n=2;n<10000000;n++)
{
p=1;
G=0;
X=0;
V=(int)ceil(pow(((float)n-X),0.5));
//printf("\nV=%d\n",V);
for(i=2;i<=V+1;i++)
{
if((n%i)==0) p=0;
else if(G==0)
{
X=i;
G=1;
V=(int)ceil(pow(((float)n-X),0.5));
}
S++;
}

if(p==0)
;//printf("NOK");
else
fprintf(f,"%Ld\n",n);
}
printf("\nAm terminat...");
fprintf(f,"S=%Lf",S);
fclose(f);
getch();

}

Farey generation using the Denominator Recurrence

/******************************************************************
//
/* Generation of FAREY Sequence using the Denominator recurrence */
//
/******************************************************************
/*
f(2*n+1) = f{n) f(2*n+2} = f(n) + f(n+1),

n=0, 1, 2, ...
*/
/******************************************************************


#include "stdio.h"
#include "math.h"
#include "conio.h"

typedef unsigned __int64 UL_int;

UL_int i, n;
UL_int nom, denom;

UL_int MAX = 0;

UL_int SIR[(0x7FFFFFFF>>4)];

FILE * pf;

UL_int f(UL_int n)
{
    if(SIR[n]!=0) return SIR[n];
    else
    {


if ( (n==1) || (n==2) ) {SIR[n] = 1; return 1;}
else if((n%2)==1) {SIR[n] = (f((UL_int)floor(((n-1)/2.0))) + f((UL_int)floor((n+1)/2.0))); return SIR[n];}
else {SIR[n] = f((UL_int)floor(n/ 2.0)); return SIR[n];}
    }
}

int main(void) {
//clrscr();

pf=fopen("fractions.txt","w+");

//printf ("\nREC: -> NOM(n) / DENOM(n):\n");
for(i=1;i<8388608 16="16" 2="2" bits="bits" for="for" i="i" numbers="numbers" on="on" p="p">{
nom = f (i) ;
denom = f (i+1) ;

if(MAXif(MAX
//printf ("\ni=%llu %llu / %llu",i,nom,denom) ;
//fprintf (pf,"\ni=%llu %llu / %llu",i,nom,denom) ;
}

printf ("\nMAX= %llu", MAX) ;

getch () ;

fclose(pf);

return 1;
}

Reference:
http://www.cut-the-knot.org/blue/Fusc.shtml

Chord - Tangent method to solve equations

//****************************************************************************
# include ;
# include ;
# include ;

# define EPS 1.e-6

double tc,tp,sc,sp,xp,xc,a,b,aux;
long int NrIt;


//****************************************************************************
double f(double x)
{
return x*x*x*x-x*x-2;
}

double df(double x)
{
return 4*x*x*x-2*x;;
}

//****************************************************************************


void main(void)
{
clrscr();
cout<<"Please enter the limits of the interval :\n";
cin>>a>>b;

NrIt=0;

if(f(a)*df(a)>0)
{
tp=a;
sp=b;
}
else if(f(a)*df(a)<0)
{
tp=b;
sp=a;
}

aux=0;
xp=xc=(sp+tp)/2.0;

while((f(xc)>=EPS)/*&&(fabs(xc-aux)>=EPS)*/)
{
NrIt++;

tc=tp-f(tp)/df(tp);
sc=(sp*f(tp)-tp*f(sp))/(f(tp)-f(sp));

xc=(sc+tc)/2.0;
aux=xp;
xp=xc;
tp=tc;
sp=sc;

cout<<"\nX"<
}

cout<<"\nThe solution is: "<cout<<"\n\nNo. of iterations= "<
getche();
}

Farey Chord - Tangent method to solve equations

If the limits of the interval are very big than the Chord-Tangent method using the Farey addition produce in half of steps better approximations:


//****************************************************************************
# include ;
# include ;
# include ;

# define EPS 1.e-6

double tc,tp,sc,sp,xp,xc,a,b,aux,mij2;
double v_N[3],v_D[3];//[0] pt a, [1] pt mij, [2] pt b
long int NrIt=0;


//****************************************************************************
double f(double x)
{
return x*x*x*x-x*x-2;
}

double df(double x)
{
return 4*x*x*x-2*x;;
}

//****************************************************************************


void main(void)
{
clrscr();
cout<<"Please enter the limits of the interval :\n";
cin>>a>>b;

NrIt=0;

if(f(a)*df(a)>0)
{
tp=a;
sp=b;
v_N[0]=b;//sp
v_D[0]=1;
v_N[2]=a;//tp
v_D[2]=1;
}
else if(f(a)*df(a)<0)
{
tp=b;
sp=a;
v_N[0]=a;//sp
v_D[0]=1;
v_N[2]=b;//tp
v_D[2]=1;
}

//double mij;

v_N[1]=v_N[0]+v_N[2];
v_D[1]=v_D[0]+v_D[2];
xp=xc=v_N[1]/v_D[1];//pt xp, xc

//mij=
aux=0;
mij2=(a+b)/2;

while(
(fabs(f(xc))>=EPS)
&&(fabs(f(xc)/xc)>=EPS)
&&(fabs(xc-aux)>=EPS)
)
{
NrIt++;

tc=tp-f(tp)/df(tp);
sc=(sp*f(tp)-tp*f(sp))/(f(tp)-f(sp));

xc=(v_N[0]+v_N[2])/(v_D[0]+v_D[2]);

if(f(sc)*f(xc)<0)
{
if(mij2>xc)
{
if(f(sc)*f(xc)<0)
{
v_N[2]=v_N[1];
v_D[2]=v_D[1];
tc=v_N[2]/v_D[2];
}
}
else //mij2<=xc
{
if(f(sc)*f(mij2)<0)
{
v_N[2]=mij2;
v_D[2]=1;
tc=v_N[2]/v_D[2];
}
else if(f(mij2)*f(xc)<0)
{
v_N[0]=mij2;
v_D[0]=1;
sc=v_N[2]/v_D[2];
v_N[2]=v_N[1];
v_D[2]=v_D[1];
tc=v_N[2]/v_D[2];
}
}
}

else if(f(xc)*f(tc)<0)
{
if(mij2 {
if(f(mij2)*f(xc)<0)
{
v_N[0]=mij2;
v_D[0]=1;
sc=v_N[0]/v_D[0];
}
}
else //mij2>=xc
{
if(f(xc)*f(mij2)<0)
{
v_N[0]=v_N[1];
v_D[0]=v_D[1];
sc=v_N[0]/v_D[0];
v_N[2]=mij2;
v_D[2]=1;
tc=v_N[0]/v_D[0];
}
else if(f(mij2)*f(tc)<0)
{
v_N[0]=mij2;
v_D[0]=1;
sc=v_N[0]/v_D[0];
}
}
}



v_N[1]=v_N[0]+v_N[2];
v_D[1]=v_D[0]+v_D[2];
xc=(v_N[0]+v_N[2])/(v_D[0]+v_D[2]);

mij2=(sc+tc)/2;

aux=xp;
xp=xc;
tp=tc;
sp=sc;

cout<<"\nX"<}

cout<<"\nThe solution is: "<cout<<"\n\nNo. of iterations= "<
getche();
}

Program Halley - Farey to solve equations

//****************************************************************************
// Program Halley, Bisectie, BisFarey, Hall_Bis, Hall_BissFarey
//****************************************************************************
# include ;
# include ;
# include ;
# include ;
# include ;

# define epx 1e-3 //epx=eps pt F
# define EPS 0.000001
# define max 1500 //nr max de iteratii

//****************************************************************************

double F(double x) { return x*x*x*x-x*x-2; }
double dF(double x) { return 4*x*x*x-2*x; }
double ddF(double x) { return 12*x*x-2; }

//****************************************************************************

double a,aux,a_aux,b_aux,b,mij,epx_aux,val,v_N[3],v_D[3];
//[0] pt a, [1] pt mij, [2] pt b
int c,opt; //contor

double x0,xp,xc;
long int NrIt;


//****************************************************************************
void main()
{
clrscr();
opt=0;

while(opt!=6)
{
clrscr();
cout<<"Please choose an option :\n";
cout<<"\t1. Method bisection\n";
cout<<"\t2. Metoda bisection using the FAREY addition\n";

cout<<"\t3. Method Halley\n";
cout<<"\t4. Method Halley_Bisection\n";
cout<<"\t5. Method Halley_BisectionFarey\n";

textcolor(YELLOW);
cprintf(" 6. Enter the interval");
textcolor(15);
cout<<"\n\t7. END THE PROGRAM\n";
cin>>opt;
if(opt!=6) cout<<"\a";
}

while(opt<=6)
{
c=1;
if(opt==6)
{
cout<<"Enter the limit of the interval, a,b :\n";
cin>>a>>b;
a_aux=a;
b_aux=b;
}

if(opt<=2)
{
aux=0;
a=a_aux;
b=b_aux;
mij=(a+b)/2;

//initializari
if(opt==1) mij=(a+b)/2;
else if(opt==2)
{
v_N[0]=a;
v_D[0]=1;
v_N[2]=b;
v_D[2]=1;
v_N[1]=v_N[0]+v_N[2];
v_D[1]=v_D[0]+v_D[2];
mij=v_N[1]/+v_D[1];
}

c=0;
epx_aux=epx+1;
val=fabs(F(mij));

while((c=epx)||(fabs(aux-mij)>epx)||(fabs(val/mij)>epx))
{
epx_aux=fabs(F(mij)/mij);

if(F(a)*F(mij)<0)
{
if(opt==1) b=mij;
else
{
v_N[2]=v_N[1];
v_D[2]=v_D[1];
b=v_N[2]/v_D[2];
}
}
else if(F(mij)*F(b)<0)
{
if(opt==1) a=mij;
else
{
v_N[0]=v_N[1];
v_D[0]=v_D[1];
a=v_N[0]/v_D[0];
}
}

aux=mij;
if(opt==1) mij=(a+b)/2;
else
{
v_N[1]=v_N[0]+v_N[2];
v_D[1]=v_D[0]+v_D[2];
mij=(v_N[0]+v_N[2])/(v_D[0]+v_D[2]);
}

c++;

cout.precision(30);
cout<< endl <<" Epx="<< epx_aux<<" ; nr it="<< c;
val=fabs(F(mij));
}
cout.precision(10);
cout<<"\n\nAm iesit pe :";

if(c>=max) cout<<"\n No. of iterations ="<< c;
if(fabs(F(mij))<=epx) cout<<"\n Criterion with epx ="<< fabs(F(mij)/(mij));

cout<<"\n";
cout<<"\nNo. of iterations PRG ="<< max;

cout<< endl<< endl<<"Solution of the equation is, x= "<< mij;

getche();
}

if(opt==3)
{
clrscr();
cout<<"Enter the start value: ";
cin>>x0;

NrIt=0;

xp=x0;
xc=xp;
aux=0;
while(fabs((F(xc))>=EPS)||(fabs(aux-xc)>=EPS))
{
NrIt++;
xc=xp-(2*F(xp)*dF(xp))/(2*dF(xp)*dF(xp)-F(xp)*ddF(xp));
aux=xp;
xp=xc;
}

cout<<"\nThe solution is: "<< xc<<" with the precision "<< EPS;
cout<<"\n\nNo. of iterations= "<< NrIt;

getche();

}

//****************************************************************************

if(opt==4)
{
clrscr();
cout<<"Enter the start value : ";
cin>>x0;

NrIt=0;

double _xa,_xb,_mij,_aux,DIF;

xp=x0;
xc=xp;
aux=0;

_xa=a;
_xb=b;
_mij=(a+b)/2;
_aux=0;
while(fabs((F(xc))>=EPS)||(fabs(aux-xc)>=EPS)||(fabs(_aux-xc)>=EPS))
{
NrIt++;

//
if(F(_xa)*F(_mij)<0)_xb=_mij;
else if(F(_xb)*F(_mij)<0) _xa=_mij;
_mij=(_xa+_xb)/2;
DIF=_mij-_aux;
_aux=_mij;
//

xc=xp-(2*F(xp))/(2*dF(xp)-ddF(xp)*DIF);
aux=xp;
xp=xc;
}

cout<<"\nThe solution is: "<< xc<<" with the precision "<< EPS;
cout<<"\n\nNo of iteartions= "<< NrIt;

getche();

}

//****************************************************************************

if(opt==5)
{
clrscr();
cout<<"Enter the start value : ";
cin>>x0;

NrIt=0;

double _xa,_xb,_mij,_aux,DIF;

xp=x0;
xc=xp;
aux=0;

_xa=a;
_xb=b;
//_mij=(a+b)/2;
//_aux=mij;

v_N[0]=a;
v_D[0]=1;
v_N[2]=b;
v_D[2]=1;
v_N[1]=v_N[0]+v_N[2];
v_D[1]=v_D[0]+v_D[2];
_mij=v_N[1]/+v_D[1];
_aux=0;

while(fabs((F(xc))>=EPS)||(fabs(aux-xc)>=EPS)||(fabs(_aux-xc)>=EPS))
{
NrIt++;

//
if(F(_xa)*F(_mij)<0)
{
v_N[2]=v_N[1];
v_D[2]=v_D[1];
_xb=v_N[2]/v_D[2];
}
else if(F(_xb)*F(_mij)<0)
{
v_N[0]=v_N[1];
v_D[0]=v_D[1];
_xa=v_N[0]/v_D[0];
}

v_N[1]=v_N[0]+v_N[2];
v_D[1]=v_D[0]+v_D[2];
_mij=(v_N[0]+v_N[2])/(v_D[0]+v_D[2]);

DIF=_mij-_aux;
_aux=_mij;
//

xc=xp-(2*F(xp))/(2*dF(xp)-ddF(xp)*DIF);
aux=xp;
xp=xc;
}

cout<<"\nThw solution is: "<< xc<<" with the precision "<< EPS;
cout<<"\n\nNo. of iterations= "<< NrIt;

getche();

}


//****************************************************************************
clrscr();

cout<<"Please choose an option :\n";
cout<<"\t1. Method bisection\n";
cout<<"\t2. Metoda bisection using the FAREY addition\n";

cout<<"\t3. Method Halley\n";
cout<<"\t4. Method Halley_Bisection\n";
cout<<"\t5. Method Halley_BisectionFarey\n";

printf("\t6. Enter the interval");
cout<<"\n\t7. END THE PROGRAM\n";

cin>>opt;
}

}
//****************************************************************************

Program Bisection method - Farey addition

# include ;
# include ;
# include ;
# include ;
# include ;
# include ;

# define epx 1e-17 //epx=eps pt F
# define max 1500 //nr max de iteratii

long unsigned int cmmdc(long unsigned int a,long unsigned int b)
{
long unsigned int aux;
while(fmod(a,b))
{
aux=a;
a=b;
b=fmod(aux,a);
}
return b;
}

//-------------------------------------
long double F(long double x) { return x*x-3; }
//-------------------------------------

long double a,a_aux,b_aux,b,dc,mij,epx_aux,val;
long unsigned int ia,ia_aux,ib_aux,ib,v_N[3],v_D[3];
//[0] pt a, [1] pt mij, [2] pt b

int c,opt,pa,pb,pmax; //contor,optiune,puteri

void main()
{
clrscr();
opt=0;

while(opt!=3)
{
clrscr();
cout<<"Alegeti o optiune :\n";
cout<<"\t1. Metoda bisectiei\n";
cout<<"\t2. Metoda bisectiei utilizand insumarea FAREY\n";
cout<<"\t3. Introducere interval\n";
cout<<"\t4. TERMINARE PROGRAM\n";
cin>>opt;
if(opt!=3) cout<<"\a";
}

while(opt<=3)
{
c=1;
if(opt==3)
{
cout<<"Introduceti a,b :\n";

//Conversie la tipul: unsigned long int
char * sa,* sb;
int poz;
gets(sa);
gets(sb);
a=a_aux=_atold(sa);
b=b_aux=_atold(sb);

int i;

//pentru a --------------
poz=0; //-daca nr e intreg
for(i=0;i < strlen(sa);i++)
if(sa[i]=='.') poz=i;
pa=strlen(sa)-poz-1;//exponentul lui a; (puterea lui a)

//pentru b --------------
poz=0; //-daca nr e intreg
for(i=0;i if(sb[i]=='.') poz=i;
pb=strlen(sb)-poz-1;//exponentul lui b; (puterea lui b)

//Calcul putere maxima
if(pa>pb) pmax=pa;
else pmax=pb;
//

ib=b*pow10(pmax);
ia=a*pow10(pmax);

//Afisare de control
cout.precision(12);
cout<<"Am citit datele:\n"<< a <<" "<< b;
cout<<"\nDatele convertite sunt:\n"<< ia <<" "<< ib;
getche();
}

if(opt<=2)
{
a=a_aux;
b=b_aux;

//Initializari
if(opt==1) mij=(a+b)/2;
else if(opt==2)
{
v_N[0]=ia;
v_D[0]=pow10(pmax);
//----------------------
dc=cmmdc(v_N[0],v_D[0]);
cout<<"\n v_*[0]= "< v_N[0]=v_N[0]/dc;
v_D[0]=v_D[0]/dc;
cout<<"Dupa siplif cu "< cout<<"v_*[0]= "< getche();
//----------------------
v_N[2]=ib;
v_D[2]=pow10(pmax);
//----------------------
dc=cmmdc(v_N[2],v_D[2]);
cout<<"\n v_*[2]= "< v_N[2]=v_N[2]/dc;
v_D[2]=v_D[2]/dc;
cout<<"Dupa siplif cu "< cout<<"v_*[2]= "< getche();
//----------------------
v_N[1]=v_N[0]+v_N[2];
v_D[1]=v_D[0]+v_D[2];
//----------------------
dc=cmmdc(v_N[1],v_D[1]);
cout<<"\n v_*[1]= "< v_N[1]=v_N[1]/dc;
v_D[1]=v_D[1]/dc;
cout<<"Dupa siplif cu "< cout<<"v_*[1]= "< getche();
//----------------------
mij=(long double)(v_N[1])/(long double)(v_D[1]);
}

c=0;
epx_aux=epx+1;
val=fabs(F(mij));

while((c < max)&&(val >= epx))
{
epx_aux=fabs(F(mij)/mij);

if(F(a)*F(mij)< 0)
{
if(opt==1) b=mij;
else
{
v_N[2]=v_N[1];
v_D[2]=v_D[1];
//----------------------
dc=cmmdc(v_N[2],v_D[2]);
cout<<"\n v_*[2]= "< v_N[2]=v_N[2]/dc;
v_D[2]=v_D[2]/dc;
cout<<"Dupa siplif cu "< cout<<"v_*[2]= "< getche();
//----------------------
b=(long double)v_N[2]/(long double)v_D[2];
}
}
else if((F(mij)*F(b)< 0)&&(F(a)*F(mij) > 0))
{
if(opt==1) a=mij;
else
{
v_N[0]=v_N[1];
v_D[0]=v_D[1];
//----------------------
dc=cmmdc(v_N[0],v_D[0]);
cout<<"\n v_*[0]= "< v_N[0]=v_N[0]/dc;
v_D[0]=v_D[0]/dc;
cout<<"Dupa siplif cu "< cout<<"v_*[0]= "< getche();
//----------------------
a=(long double)v_N[0]/(long double)v_D[0];
}
}

if(opt==1) mij=(a+b)/2;
else
{
v_N[1]=v_N[0]+v_N[2];
v_D[1]=v_D[0]+v_D[2];
//----------------------
dc=cmmdc(v_N[1],v_D[1]);
cout<<"\n v_*[1]= "< v_N[1]=v_N[1]/dc;
v_D[1]=v_D[1]/dc;
cout<<"Dupa siplif cu "< cout<<"v_*[1]= "< //----------------------
mij=(long double)(v_N[1])/(long double)(v_D[1]);
}

c++;

cout.precision(12);
cout<
if(c>=max) cout<<"\n Nr. de it. ="< if(fabs(F(mij))<=epx) cout<<"\n Criteriu cu epx ="<
cout<<"\n";
cout<<"\nNr. de it. PRG ="<
cout< cout<<"\t1. Metoda bisectiei\n";
cout<<"\t2. Metoda bisectiei utilizand insumarea FAREY\n";
cout<<"\t3. Modificare interval\n";
cout<<"\t4. TERMINARE PROGRAM\n";
cin>>opt;

}

}