/*
*  Fibonacci
*
*  Author	: Rosario Turco
*  Date		: 31/12/2009 
*  Revision	: 1.0           
*/

/* abbreviazione */
Fib(n)= local(); { 
 return(fibonacci(n););
}

/*
*  
*  0 - Nok
*  1 - Ok   
*
*/
isFib(n) = local(val=0, status=0); { 
 val = 5 * n^2 + 4;
 if( frac(n) == 0.0 & issquare(val) == 1, status = 1; );
 return(status);  
}

/*
*  
*  index of Fibonacci's number   
*  -1 Nok
*/

indFib(n)= local(ind=-1); { 
 if( isFib(n) == 0, return(ind);); 
 if ( n < 0, ind = 0; return(ind););

 ind = ceil(log (sqrt (5) * abs(n)) / log (Phi()));
 return(ind);
}

/*
*  
*  Fibonacci's number complex   
*  
*/
CFib(z) = local(phi=0,F=0); { 
 phi = Phi();
 F = phi^z – cos(z*Pi);
 return(F / (phi^z) / sqrt(5););
}


/* 
*
*  Golden ratio 
*
*/
Phi()=local(); { 

 return((1 + sqrt (5)) / 2;);

}

KF(num)=local(F=0,P=1,i=1); { 

 F = sum(i=1, num, 1.0/fibonacci(i));
 num=num+1;

 while( i != num,
        a = fibonacci(i);
        b = 1.0/a;
        if(isprime(a) == 1,
           P = P * (1-b);
        );
        i=i+1;   
  );

  return(F * P);
}


Fibor(num=100)=local(P=1,i=1, a=1); { 

 num=num+1;

 while( i != num,
        a = fibonacci(i);
        P = P * a;
        i=i+1;   
  );

  return(P);
}



/*
*   Keith's numbers Research (Repfigit and RevRepfigit)
*   R. Turco
*   rep=1 if I search Repfigit 
*   rep=0 if I search RevRepfigit
*
*/

KNRA(start=10, end=10000000, rep=1)=local();{ 
 
 for(i=start,end,
      KNR(i,rep);   
 );

}

KNR(val,rep)=local(value=0,p=0,s="",str="", out="", m=0, j=0, s0=0, s1=0, s2=0, trovato=0); { 

  
     s = concat(s,val);       
     m = length(s);

     MyW=Vecsmall(s);
     MyW2=Vecsmall(s);
     v=vector(m);
     v2=vector(m);

     /* Reverse Vector */ 

     MyW2=RevVect(s,m); 
          

     /* I calculate the sum */   

     for( j=1,m,
          v[j] = MyW[j]-48;          
          v2[j] = MyW2[j]-48;
          out = concat( out, v[j] );
          out = concat( out, "," );
          str = concat(str, v2[j]);
          s1 = s1 + v[j];
     );
     out = concat( out, s1 );
     out = concat( out, "," );

     v[1]=s1; 

     /* 
         I calculate the first value of the vector.
         This is important in the revRepfigit 
     */

     value=GetValStr(s,m);

     /* I get the value of the reverse vector */

     s2=GetValStr(str,m);
     
     p=1;

     /* I search the Repfigit */
     if( rep == 1,   

         while( trovato == 0 & s1<val+1, 
            p++;
            if( p>m, p=1;);
            s1 = 0;

            for( j=1,m,
              s1 = s1 + v[j]; 
            );
 
            out = concat( out, s1 );
            if( s1 == val, print("Repfigit trovato : ", val, " len : ", m); trovato = 1; print(out););
            if( s1 != val, out = concat( out, "," ););
            v[p] = s1;
             
        );
     );

     /* revRepfigit */

     if( rep == 0 & (10 < s2) & (value%10)!=0,   
 
         while( trovato == 0 & s1<s2+1, 
            p++;
            if( p>m, p=1;);
            s1 = 0;

            for( j=1,m,
              s1 = s1 + v[j]; 
            );
 
            out = concat( out, s1 );
            if( s2 == s1, print("revRepfigit trovato : ", val, " len : ", m); trovato = 1; print(out);); 
            if( s1 != val, out = concat( out, "," ););
            v[p] = s1;
             
        );
    );

 
}



GetValStr(str,len)=local(val=0,weight=0);{ 

 MyStr=Vecsmall(str);
 weight=len-1; 
 for(k=1,len,
         val = (MyStr[k]-48)*(10^(weight)) + val;
         weight--; 
 );
 return(val);
}

RevVect(str,len)=local(k=0,j=0);{ 

     MyWA=Vecsmall(str);
     MyWB=Vecsmall(str);
 
     k=len;
     j=1;
     while( k>0,
            MyWB[k]=MyWA[j];
            k--;
            j++;
     );

     return(MyWB);
}