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.

313 lines
8.5 KiB

#!/usr/local/bin/bc -l
### CF.BC - Continued fraction experimentation using array pass by reference
# Best usage for the output functions in this library:
# .=output_function(params)+newline()
# This suppresses the output and also appends a newline which is not added
# by default.
## Initialisation / Tools / Workhorses ##
# Borrowed from funcs.bc
define int(x) { auto os;os=scale;scale=0;x/=1;scale=os;return(x) }
# Borrowed from output_formatting.bc
define newline() { print "\n" }
# sanity check for the below
define check_cf_max_() {
auto maxarray;
maxarray = 4^8-4;cf_max=int(cf_max)
if(1>cf_max||cf_max>maxarray)cf_max=maxarray
return 0;
}
# global var; set to halt output and calculations at a certain point
cf_max=-1; # -1 to flag to other libs that cf.bc exists before value is set
# Workhorse function to prepare for and tidy up after CF generation
define cf_tidy_(){
# POSIX scope; expects vars cf__[], max, p and i
# Tidy up the end of the CF
# assumes the last element of large CFs to be invalid due to rounding
# and deletes that term.
# for apparently infinite (rather than just long) CFs, uses a bc trick
# to signify special behaviour to the print functions.
# for apparently finite CFs, checks whether the CF ends ..., n, 1]
# which can be simplified further to , n+1]
auto j,diff,rl,maxrl,bestdiff
cf__[i]=0
if(p>max){
.=i--;
if(p<max*A^5){ # assume infinite
cf__[i]=0.0; # bc can tell between 0 and 0.0
# Identify repeating CFs and store extra information
# that is used by the print functions.
max=i;mrl=0;
p=max-1;if((i=p-1)>0)for(i=i;i;i--){
if(cf__[i]==cf__[p]){
diff=p-i;rl=0
for(j=p-1;j>diff;j--){if(cf__[j-diff]!=cf__[j])break;.=rl++}
if(j<=i)rl+=diff
if(rl>mrl){mrl=rl;bestdiff=diff}
}
}
if(3*mrl>=p){cf__[++max]=p-mrl-1;cf__[++max]=bestdiff;cf__[0]+=0.0}
return 0;
}
cf__[i]=0; # bc can differentiate between 0 and 0.0
}
if(i<2)return 0;
.=i--;if(abs(cf__[i])==1){cf__[i-1]+=cf__[i];cf__[i]=0}
return 0;
}
# Workhorse function for cf_new and cfn_new
define cf_new_(near) {
# POSIX scope; expects array *cf__[] and var x
auto os, i, p, max, h, n, temp;
.=check_cf_max_()
os=scale;if(scale<scale(x))scale=scale(x)
if(scale<6+near){
print "cf";if(near)print"n";print"_new: scale is ",scale,". Poor results likely.\n"
}
max=A^int(scale/2);p=1
if(near)h=1/2;
for(i=0;i<cf_max&&p<max;i++) {
if(near){
n=int(temp=x+h);if(n>temp).=n-- # n=floor(x+.5)
} else {
n=int(x)
}
p*=1+abs(cf__[i]=n);
x-=cf__[i]
if(x==0||p==0){.=i++;break}
x=1/x
}
scale=os
return cf_tidy_();
}
## Making/unmaking continued fractions and ordinary fractions ##
# Create a continued fraction representation of x in pbr array cf__[]
define cf_new(*cf__[],x) {
return x+cf_new_(0);
}
# Create a continued fraction representation of x in pbr array cf__[]
# using signed terms to guarantee largest magnitude following term
define cfn_new(*cf__[],x) {
return x+cf_new_(1);
}
# Copy a continued fraction into pbr array cfnew__[] from pbv array cf__[]
define cf_copy(*cfnew__[],cf__[]){
auto e,i;
for(i=0;i<=cf_max&&e=cf__[i];i++)cfnew__[i]=e
if(i<=cf_max)cfnew__[i]=cf__[i] # might be 0.0
e=cf__[0];if(int(e)==e&&scale(e)){
# copy extra info
if(++i<=cf_max)cfnew__[i]=cf__[i]
if(++i<=cf_max)cfnew__[i]=cf__[i]
}
return (i>cf_max);
}
# Convert pbv array cf__[] into a number
define cf_value(cf__[]) {
auto n, d, temp, i;
.=check_cf_max_();
if(cf__[1]==0)return cf__[0];
for(i=1;i<cf_max&&cf__[i];i++){}
n=cf__[--i];d=1
for(i--;i>=0;i--){temp=d;d=n;n=temp+cf__[i]*d}
return(n/d);
}
# Convert pbv array cf__[] into a new cf of the other type
# . with respect to nearness.
define cf_toggle(*cf__[],cf2__[]){
auto os,p,i,x,zero,sign,near;
sign=1;if(cf2__[0]<0)sign=-1
near=0
for(i=1;x=cf2__[i];i++){
p*=1+abs(x)
if(x*sign<0)near=1
}
zero=x;
os=scale;scale=0
for(i=2;p;i++)p/=A
if(i<os)i=os
scale=i
x=cf_value(cf2__[])
.=cf_new_(1-near)
scale=os
for(i=0;cf__[i];i++){}
cf__[i]=zero;
return 0
}
# Return the value of the specified convergent of a CF
define cf_get_convergent(cf__[],c){
auto ocm,v;
if(c==0)return cf_value(cf__[]);
.=check_cf_max_();ocm=cf_max
cf_max=c;v=cf_value(cf[])
cf_max=ocm;return v;
}
# Return the value of the specified convergent of the CF of x
define get_convergent(x,c){
auto cf[];
.=cf_new(cf[],x)
return cf_get_convergent(cf[],c)
}
# Create denominator, numerator and optional intpart in
# . first three elements of pbr array f__[]
# . from CF in pbv array cf[]()
# NB: returned elements are in reverse of the expected order!
define frac_from_cf(*f__[],cf__[],improper) {
auto n,d,i,temp;
.=check_cf_max_();
improper=!!improper;
if(cf__[0]==(i=int(cf__[0])))cf__[0]=i
if(cf__[1]==0){f__[0]=1;f__[1]=0;return f__[2]=cf__[0]}
for(i=1;i<cf_max&&cf__[i];i++){}
n=cf__[--i];d=1
for(i--;i>=!improper;i--){
temp=n;n=d;d=temp # reciprocal = swap numerator and denominator
n+=cf__[i]*d
}
temp=0
if(!improper){temp=n;n=d;d=temp}#correct for having stopped early
if(d<0){n=-n;d=-d} # denominator always +ve
if(!improper&&cf__[0]!=0){
temp=cf__[0]
}
f__[0]=d;f__[1]=n;f__[2]=temp
return temp+n/d;
}
# Upscale an allegedly rational number to the current scale
define upscale_rational(x) {
auto os,f[],cf[];
if(scale<=scale(x))return x
os=scale;scale=scale(x)
.=cf_new(cf[],x); # Sneaky trick to upscale (literally)
.=frac_from_cf(f[],cf[],1) # any rational value of x
scale=os
x=f[1]/f[0] # x is now (hopefully) double accuracy!
return x;
}
## Output ##
# Set to 1 to truncate quadratic surd output before repeat occurs
cf_shortsurd_=0
# Output pbv array cf__[] formatted as a CF
define cf_print(cf__[]) {
auto i,sign,surd,sli,sri;
.=check_cf_max_();
sign=1;if(cf__[0]<0||(cf__[0]==0&&cf__[1]<0)){sign=-1;print "-"}
# Check for surd flag and find lead in and repeat period
surd=0;sri=int(i=cf__[0]);if(scale(i)&&sri==i){
for(i=1;cf__[i];i++){}
sli=cf__[++i]
surd=cf__[++i]
}
print "[", sign*sri;
if(cf__[1]==0){print "]";return 0}
print "; ";
if(surd)if(sli){sri=sli-1}else{sri=surd-1}
for(i=2;i<=cf_max&&cf__[i];i++){
print sign*cf__[i-1]
if(!surd||sri){
print ", ";if(surd).=sri--
} else {
if(sli){
sli=0
}else if(cf_shortsurd_){print ";...]";return 0}
print "; ";sri=surd-1
}
}
print sign*cf__[i-1];
# detect a 0.0 which signifies that the cf has been determined
# to be infinite (not 100% accurate, but good)
if(scale(cf__[i])){
if(!surd||sri){print ","}else{print ";"}
print "..."
}
print "]";
return (i>cf_max);
}
# Print a number as a continued fraction
define print_as_cf(x) {
auto cf[];
.=cf_new(cf[],x)+cf_print(cf[])
return x;
}
# Print a number as a signed continued fraction
define print_as_cfn(x) {
auto cf[];
.=cfn_new(cf[],x)+cf_print(cf[])
return x;
}
# Output pbv array cf__[] as a fraction or combination of int and fraction
# . the 'improper' parameter makes the choice.
# . . set to non-zero for top-heavy / improper fractions
define cf_print_frac(cf__[],improper) {
auto f[],v;
v=frac_from_cf(f[],cf__[],improper)
if(f[1]==0){print f[2];return f[2]}
if(!improper&&cf__[0]!=0){
print cf__[0]
if((f[1]<0)==(f[0]<0)){print"+"} # if n and d have same sign...
}
print f[1],"/",f[0] # n/d
return v;
}
# Print x as a fraction or combination of int and fraction
# . the 'improper' parameter makes the choice.
# . . set to non-zero for top-heavy / improper fractions
define print_frac(x,improper) {
auto cf[]
.=cf_new(cf[],x)+cf_print_frac(cf[],improper)
return x;
}
# Print x as a fraction or combination of int and fraction
# . the 'improper' parameter makes the choice.
# . . set to non-zero for top-heavy / improper fractions
# . This alternative function rounds to the nearest integer
# . . for the integer part and then shows the addition
# . . or _subtraction_ of the fraction.
# . . e.g. 1+2/3 is shown as 2-1/3
define print_frac2(x,improper) {
auto cf[];
.=cfn_new(cf[],x)+cf_print_frac(cf[],improper)
return x;
}
define cf_print_convergent(cf__[],c){
auto ocm,n,v;
if(c==0)return cf_print_frac(cf__[],0);
if(c<0){n=1;c=-c}else{n=c}
.=check_cf_max_();ocm=cf_max;
for(n=n;n<=c;n++){cf_max=n;v=cf_print_frac(cf__[],1);if(n!=c)print "\n"}
cf_max=ocm;return v;
}
define print_convergent(x,c){
auto cf[],v;
v=cf_new(cf[],x)
v=cf_print_convergent(cf[],c)
return v;
}