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.

92 lines
2.3 KiB

#!/usr/local/bin/bc -l funcs.bc cf.bc
### CF_SYLVESTER.BC - Sylvester expansion experimentation using array pass by reference
# Workhorse: Create an Sylvester expansion of x in pbr array sy__[]
define sylvester_new_(mode) {
# expects *sy__[] and x to be declared elsewhere
auto i,max,q,iq,p,h
# error checking
.=check_cf_max_()
#p=(abs(x)>=1);
q=(scale<6);
#if(p||q){print "sylvester";if(mode)print "2";print "_new: "}
#if(p){print "Can't work with integer part.\n"; x-=int(x)}
if(q){print "scale is ",scale,". Poor results likely.\n"}
max=A^int(scale/2)
i=0;p=1;h=3/2
for(i=0;x&&p<max&&i<cf_max;i++){
q=1/abs(x)
if(!mode){
q=ceil(q) # proper sylvester expansion
} else {
q=floor(q+h) # secondary sylvester expansion
}
p*=q
sy__[i]=q*=sgn(x)
x-=1/q
}
if(p>=max){
if(!mode)if(!--i)i=1
if(abs(sy__[i])>max){
sy__[i]=0;.=i--
while(iq=int(q=(sy__[i]*sy__[i-1])/(sy__[i]+sy__[i-1]))==q){
sy__[i]=0;sy__[--i]=iq
}
}else{
sy__[i]=0.0
}
}else{
sy__[i]=0
}
return 0;
}
#echo 'scale=250;sylvester_new(a[],7/15);.=sylvester_print(a[])+newline();sylvester_value(a[])' | bc -l funcs.bc *.bc
# Create Sylvester expansion of x in pbr array sy__[]
# . all terms same sign as x
define sylvester_new(*sy__[],x) {
return x+sylvester_new_(0);
}
# Create secondary Sylvester expansion of x in pbr array sy__[]
# . first term same sign as x, all other terms opposite sign
# . terms in implied Egyptian fraction sum have alternating signs
define sylvester2_new(*sy__[],x) {
return x+sylvester_new_(1);
}
# Output pbv array sy__[] formatted as an Sylvester expansion
define sylvester_print(sy__[]){
auto i;
.=check_cf_max_();
print "{";
if(sy__[1]==0){print sy__[0],"}";return 0}
for(i=1;sy__[i];i++)print sy__[i-1],", ";
print sy__[i-1];
if(scale(sy__[i]))print ",...";
print "}";return 0;
}
# Print a number as an Sylvester expansion
define print_as_sylvester(x){
auto sy[];
.=sylvester_new(sy[],x)+sylvester_print(sy[])
return x;
}
define print_as_sylvester2(x){
auto sy[];
.=sylvester2_new(sy[],x)+sylvester_print(sy[])
return x;
}
# Turn the Sylvester expansion in pbv array sy__[] into its value
define sylvester_value(sy__[]) {
auto i,p,n,d;
.=check_cf_max_();
n=0;d=1;for(i=0;i<cf_max&&p=sy__[i];i++){n=n*p+d;d*=p}
return n/d;
}