1
0
Fork 0
Libraries for bc and dc.
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

130 lines
2.6 KiB

#!/usr/local/bin/bc -l funcs.bc
### Fibonacci.BC - Fibonacci and Lucas functions
# Requires funcs.bc
# n-th Fibonacci number over the reals
define fibonacci(n){
auto a,b,c,intn,count,fracn,s5,os
if(n==0)return 0
os=scale;scale=0;count=intn=n/1
if(n<0){
scale=os;
a=-fibonacci(-n)
if(n==intn)return a*(-1)^(-intn)
return a*c(pi()*n)
}
count+=2;
a=-1;b=1;c=0
while(--count){
c=a+b;a=b;b=c
}
scale=os;
if(n==intn)return c
fracn=n-intn
s5=sqrt(5)
a=e(fracn*l( (1+s5)/2 ))
a*=(s5*c+sqrt(5*c*c+4*(-1)^intn))/2
a=(a-c(pi()*n)/a)/s5
return a
}
# inverse of the above - cannot deal with values below 1 (except 0)
# but is accurate to within 'scale' decimal places otherwise
define inverse_fibonacci(f) {
auto a,b,c,intn,intf,fracf,s5,phinx2,eps,s5f,z5f2,lph,pi,os
if(f==0)return f
if(f<1)return 0 # avoid multivalued mess
os=scale;scale=0;intf=f/1
a=-1;b=1;c=0
for(intn=-2;c<=intf;intn++){
c=a+b;a=b;b=c
}
scale=os
if(f==a)return intn
c=a
s5=sqrt(5)
phinx2=s5*c+sqrt(5*c*c+4*(-1)^intn)
lph=l( (1+s5)/2 )
pi=pi()
s5f=s5*f
z5f2=5*f*f
a=0.5 #start guess
os+=8
for(scale=8;scale<=os;scale+=8){
b=0
eps=A^(2-scale)
while(abs(a-b)>eps){
b=a
a=s5f+sqrt(z5f2+4*c(pi*(intn+a)))
a/=phinx2
a=l(a)/lph
}
a=(a+b)/2
}
os-=8;scale=os;a/=1
return intn+a
}
# n-th Lucas number over the reals
define lucas(n){
auto a,b,c,intn,count,fracn,os
if(n==0)return 2
os=scale;scale=0;count=intn=n/1
if(n<0){
scale=os;
a=lucas(-n)
if(n==intn)return a*(-1)^(-intn)
return a*c(pi()*n)
}
count+=2;
a=3;b=-1;c=2
while(--count){
c=a+b;a=b;b=c
}
scale=os;
if(n==intn)return c
fracn=n-intn
a=e(fracn*l( (1+sqrt(5))/2 ))
a*=(c+sqrt(c*c-4*(-1)^intn))/2
a=a+c(pi()*n)/a
return a
}
# inverse of the above - inaccurate with values below 2 (except -1, 0 and 1)
# but is accurate to within 'scale' decimal places otherwise
define inverse_lucas(l) {
auto a,b,c,intn,intl,fracl,phinx2,eps,l2,lph,pi,os
if(l<-1)return -1
if(-1<=l&&l<1)return ((7-3*l)*l-A)/(2*A)
if(1<=l&&l<=2)return 2-l # avoid multivalued mess
os=scale;scale=0;intl=l/1
a=3;b=-1;c=2
for(intn=-2;c<=intl;intn++){
c=a+b;a=b;b=c
}
scale=os
if(l==a)return intn
c=a
phinx2=c+sqrt(c*c-4*(-1)^intn)
lph=l( (1+sqrt(5))/2 )
pi=pi()
l2=l*l
a=0.5 #start guess
os+=8
for(scale=8;scale<=os;scale+=8){
b=0
eps=A^(2-scale)
while(abs(a-b)>eps){
b=a
a=l+sqrt(l2-4*c(pi*(intn+a)))
a/=phinx2
a=l(a)/lph
}
}
os-=8;scale=os;a/=1
return intn+a
}