GNU bc

GNU bc is a fairly ubiquitous, useful and powerful calculator that few people seem to know much about or do anything interesting with.

It can be found on most versions of Linux and BSD, and there are at least two versions available for Windows: As a standalone executable and also available as part of Cygwin.

About this page

This page exists to raise the profile of bc, to supply pre-written code for performing precision calculations (for the most part at any rate) and investigation into various areas of mathematics, as well as being an ongoing hobby project where I show off my coding and mathematical talents... or disappointing lack thereof.

Each of the scripts / files / programs / libraries - whatever you wish to call them - contain sizeable comments explaining exactly what the file is, what it does and what each of the functions do. This is not the best way of providing documentation, but the expectation is that people should use and learn from these files, and so putting some documentation in with the code is good first step.

Below you will find slightly better documentation than in the file comments, and perhaps the occasional example, but I admit that this isn't ideal; There is the intention to set up a wiki or a forum, since there's over 250kB of bc code here; Watch this space.

A forerunner to this page was linked from the bc articles in both the English, German and Japanese Wikipedias. Welcome, Wilkommen and Irasshai to readers heading in from those places!

Files, Keywords and Functions

The files linked here contain well over 500 function definitions for GNU bc; This section should provide some sort of idea as to what kind of functions can be found in each file, beyond any hint already provided by the filename.

Before we begin...

Before downloading any of these files and to avoid any puzzled moments when reading this web page, a passing familiarity with bc is recommended. The official GNU-bc manual is well worth a read.

...some other sites of interest

Very occasionally I run across other sites with GNU bc code available for download. So far there are only three in this list, but if you know of a site, or have your own and would like to be listed here, let me know.

BC Number Theory programs
The programs found here have some overlap with funcs.bc, primes.bc and cf.bc, although the latter pale in comparison. In short, if you're doing number theory in bc, the above is the place to look if I don't have what you need.
X-bc
Provides a graphical interface to bc on X-Window supporting platforms, as well as providing various extra functions, many of which are equivalent to ones found here. Disclaimer: X-bc appears to have not been updated since 2003, and I have not used the graphical interface myself.
Cálculo numérico con bc
Marc Meléndez Schofield's Spanish language gnu-bc page has many things similar to those found on this site, as well as a few things that aren't. For example, the site serves as a tutorial to help those new to bc, and there are some innovative ideas, such as using bc to generate PPM images (supported by many graphics packages) as output, which makes this site's output_graph.bc look like a toy in comparison.

Some style notes

There are some conventions that I have tried to stick to in these files to help identify certain types of function. The main conventions are:

How to use these files

You're a busy person, and perhaps as a result you've not had time to read the official GNU-bc manual. That's understandable. All you need to do is create a directory on your computer which contains the file or files you'd like to use - put your own code in there along with them if necessary - open a command or shell window on that directory (or cd to it) then type the command

/path/to/bc -l file1.bc file2.bc file3.bc ... etc.

Of course, you'll need to use the right path to the bc executable and exchange file1, file2 and so on for actual .bc filenames. You could also put the path to the bc executable into your shell's PATH variable, removing the need to type the path, but system configuration help is outside the scope of this page.

Personally, I keep all my bc files in the same folder, and since my system is set up with the aforementioned PATH entry, all I need type is:

bc -l *.bc

...and bc launches with all my functions ready to use.

For reasons that are not entirely clear, sometimes loading files in a particular order will cause the bc interpreter to crash with an error, despite no error existing in any one file of the loading code. This then means that the above trick does not work.

To work around this, loading the files in a different order, or even loading some files more than once (!) will prevent the error from happening, for example:

bc -l [a-oq-z]*.bc p*.bc

... which loads the files in alphabetical order, skipping all files beginning with p, only to load them last, may well work when the previous example does not.

Filenames here have underscores in them so that, due to ASCII ordering, parent files load before their similarly named children when using the *.bc shortcut.

You can also take advantage of bc's standard input acceptance and perform tricks like:

echo '.=divisors_store(a[],2520)+asort0(a[])+aprint0(a[])' | bc -l primes.bc array.bc

Which, if you have those two files in particular, should spit out:

{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 15, 18, 20, 21, 24, 28, 30, \
35, 36, 40, 42, 45, 56, 60, 63, 70, 72, 84, 90, 105, 120, 126, 140, \
168, 180, 210, 252, 280, 315, 360, 420, 504, 630, 840, 1260, 2520}

The above example also goes to show some of the power of the functions which can be found on this page.

Archive of all files

For the lazy and/or curious, most of the files below can be downloaded in one ZIP archive:

This filename should update whenever a change is made to one or more of the files below.

Directory of functions and functionality

Update, February 2013: Removal of some bugs introduced in the last update as well as some minor speed improvements and modified algorithms. More new functions and some new files.

Quick navigation to main sections

array.bc
Updated for February 2013
A large number of functions for managing different array formats in bc. This file uses the undocumented pass-by-reference feature for arrays .

Function names have various suffixes or none at all for dealing with many different array formats:

aequals(a__[],b__[],count)
Boolean function; Returns whether the two arrays are equal for the number of entries specified by count.
aequals0(a__[],b__[])
Boolean function; Returns whether two zero-terminated arrays are equal
aequals2r(a__[],astart,aend,b__[],bstart,bend)
Boolean function; Returns whether the specified ranges of the two arrays are equal
aequalsb(a__[],b__[],x)
Boolean function; Returns whether two arrays, terminated by x are equal. When x is 0, this is equivalent to aequals0().
aequalsl(a__[],b__[])
Boolean function; Returns whether two arrays, whose lengths are specified as the first element of the array, are equal.
aequalsr(a__[],b__[],start,end)
Boolean function; Returns whether two arrays are equal for the given range of elements.
aappend(*a__[],aend,b__[],count)
Appends count elements of the second array to the first. Note that the index of the end of the first array must be specified.
aappend0(*a__[],b__[])
Appends the contents of the second array to the first. Both arrays are treated as zero-terminated. Zero-termination is maintained.
aappendb(*a__[],b__[],x)
Appends the contents of the second array to the first. Both arrays are expected to be terminated by whatever is specified by x.
aappendl(*a__[],b__[])
Appends the contents of the second array to the first. Both arrays should have their length as their first element
aappendr(*a__[],aend,b__[],bstart,bend)
Appends the contents of the given range of the second array to the first array. Note that the index of the end of the first array must be specified.
aappendelem(*a__[],aend,e)
Appends a single element, that is, a number to the end of an array.
aappendelem0(*a__[],e)
Appends a single element to the end of a zero-terminated array. Zero-termination is maintained.
aappendelemb(*a__[],x,e)
Appends a single element to the end of an x-terminated array.
aappendeleml(*a__[],e)
Appends a single element to the end of an array which knows its own length.
acompare(a__[],b__[],count)
acompare0(a__[],b__[])
acompare2r(a__[],astart,aend,b__[],bstart,bend)
acompareb(a__[],b__[],x)
acomparel(a__[],b__[])
acomparer(a__[],b__[],start,end)
Comparison functions; Return -1 if the first array is logically before the second in lexographic order, 0 if they are equal and 1 otherwise. Behaviour and parameters are as for the aequal...() functions.
aconv0fromr(*a__[], b__[],start,end)
Conversion function; Creates zero-terminated array a from the specified range of array b.
aconvbfromr(*a__[],x,b__[],start,end)
Conversion function; Creates array a, terminated by x, from the specified range of array b.
aconvlfromr(*a__[], b__[],start,end)
Conversion function; Creates an array whose length is the first element, a, from the specified range of array b.
aconvrfrom0(*a__[],start,end,b__[])
Conversion function; Sets the specified range of array a with the values from the zero-terminated array, b.
aconvbfrom0(*a__[],x, b__[])
Conversion function; Sets the specified range of array a with the values from the array b, which is terminated by the value specified by x.
aconvlfrom0(*a__[], b__[])
Conversion function; Creates an array whose length is the first element, a, using the values from zero-terminated array, b.
aconv0fromb(*a__[], b__[],x)
Conversion function; Creates zero-terminated array a from the values of the array b, which is terminated by the value specified by x.
aconvrfroml(*a__[],start,end,b__[])
Conversion function; Sets the specified range of array a with the values from the array b, whose length is specified by its first element.
aconv0froml(*a__[], b__[])
Conversion function; Creates zero-terminated array a from the values of the array b, whose length is specified by its first element.
aconvbfroml(*a__[],x, b__[])
Conversion function; Creates array a, which is terminated by the value in x, from the values of the array b, whose length is specified by its first element.
acopy(*a__[],b__[],count)
Overwrites the first array with count elements from the second.
acopy0(*a__[],b__[])
Turns the first array into a clone of the second, assuming both are zero-terminated arrays.
acopy2r(*a__[],astart,aend,b__[],bstart,bend)
Overwrites the specified range of the first array with the specified range of the second.
acopyb(*a__[],b__[],x)
Turns the first array into a clone of the second, assuming both are x-terminated arrays.
acopyl(*a__[],b__[])
Turns the first array into a clone of the second, assuming the second array contains the length as its first element.
acopyr(*a__[],b__[],start,end)
Overwrites those elements first array with elements of the second which are specified by the given range.
ainsertat(*a__[],acount,pos,b__[],bcount)
Experimental. Inserts elements of the second array into the first, where both array's lengths are specified in the ?count parameters, and the position within the first is given.
ainsertat0(*a__[],pos,b__[])
Experimental. Inserts the contents of the second array into the first at the given position. Both arrays are expected to be zero-terminated.
ainsertatb(*a__[],pos,b__[],x)
Experimental. As above, only the array terminator is specified by x.
ainsertatl(*a__[],pos,b__[])
Experimental. Inserts the contents of the second array into the first at the given position. Both arrays should have their length stored as their first element.
ainsertatr(*a__[],astart,aend,pos,b__[],bstart,bend)
Experimental. Inserts the specified range of the second array at the given position within the specified range of the first.
ainsertelemat(*a__[],count,pos,e)
Experimental. Inserts a single number into an array of length count at the given position.
ainsertelemat0(*a__[],pos,e)
Experimental: Inserts a single number into a zero-terminated array at the given position.
ainsertelematb(*a__[],x,pos,e)
Experimental: Inserts a single number into an x-terminated array at the given position.
ainsertelematl(*a__[],pos,e)
Experimental: Inserts a single number into an array at the given position. The array will contain its length as its first element.
ainsertelematr(*a__[],start,end,pos,e)
Experimental: Inserts a single number at the given position within the specified range of the array.
amatcharray(small__[],scount,large__[],lcount)
Returns the position, if present, of the smaller array within the larger. Array lengths are specified in the ?count parameters.
amatcharray0(small__[],large__[])
Returns the position, if present, of the smaller array within the larger. Both arrays assumed to be zero-terminated.
amatcharray2r(small__[],astart,aend,large__[],bstart,bend)
Returns the position, if present, of the given range of the smaller array within the given range of the larger.
amatcharrayb(small__[],large__[],x)
Returns the position, if present, of the smaller array within the larger. Both arrays assumed to be x-terminated.
amatcharrayl(small__[],large__[])
Returns the position, if present, of the smaller array within the larger. Both arrays assumed to have their length as their first element.
amatchelement(e,a__[],count)
Returns the position, if present, of a number within an array. Array length is specified by the count parameter.
amatchelement0(e,a__[])
Returns the position, if present, of a number within a zero-terminated array.
amatchelementb(e,a__[],x)
Returns the position, if present, of a number within an x-terminated array.
amatchelementl(e,a__[])
Returns the position, if present, of a number within an array whose length is stored as the first element.
amatchelementr(e,a__[],start,end)
Returns the position, if present, of a number within the given range of an array.
aparts3r(*a__[],astart,aend,b__[],bstart,bend,c__[],cstart,cend)
Overwrites the first array with the given ranges of itself and two further arrays, either of which may also be the original array.
aprint(*a__[],count)
Prints the first count elements of an array.
aprint0(*a__[])
Prints all elements of an array up to, but not including, the first zero (zero-terminated array).
aprintb(*a__[],x)
Prints all elements of an array up to, but not including, the first instance of x.
aprintl(*a__[])
Prints all elements (but not the length itself) of an array whose length is stored as its first element.
aprintr(*a__[],start,end)
Prints the given range of an array.
aprintu(*a__[],x)
Prints all elements of an array up to and including the first instance of x. Useful for showing the terminator on terminated arrays. x can, of course, be zero.
areverse(*a__[],n)
Reverse the first n elements of an array.
areverse0(*a__[])
Reverse a zero-terminated array.
areverseb(*a__[],x)
Reverse an x-terminated array.
areversel(*a__[])
Reverse an array whose length is stored in its first element.
areverser(*a__[],start,end)
Reverse the given range of an array.
arunlength(*v__[],*r__[], a__[], n)
Fills the first array with the values from the first n values of the third, filling the second array with an associated count (run-length) of consecutive values. Returns the length of the first two arrays (which is necessarily the same number). Note that this will be less than or equal to n. Works best on a sorted array
arunlength0(*v__[],*r__[], a__[])
Fills the first array with the values from the third array, which is zero-terminated. The second array is filled with the associated count (run-length) of consecutive values. Both of the first two arrays are then zero-terminated like the third. Returns the number of elements in these arrays (before the terminator). Works best on a sorted array
arunlengthb(*v__[],*r__[], a__[],x)
Fills the first array with the values from the third array, which is terminated by the value given in x. The second array is filled with the associated count (run-length) of consecutive values. Both of the first two arrays are then terminated like the third. Returns the number of elements in these arrays (before the terminator). Works best on a sorted array
arunlengthl(*v__[],*r__[], a__[])
Fills the first array with the values from the third array, whose length is given in its first element. The second array is filled with the associated count (run-length) of consecutive values. Both of the first two arrays are set to be in the same format as the third, having their length as their first element. Returns the number of elements in these arrays, which is the same as the aforementioned length. Works best on a sorted array
arunlengthr(*v__[],*r__[], a__[],start,end)
Fills the first array with the values from the values found in the third array between the start and end indexes inclusive. The second array is filled with the associated count (run-length) of consecutive values. Returns the length of the first two arrays. Works best on a sorted array
asanerange_(), asanerange2_()
Internal function for checking and repairing the sanity of ranges in functions that use them.
aset8(*a__[],start,a,b,c,d,e,f,g,h)
aset16(*a__[],start,a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p)
Since bc has no means of assigning to an array quickly, these functions allow the setting of 8 or 16 elements of an array in a single function call. The start parameter is the location that the first parameter will be set; The rest will then be assigned to subsequent elements.
asort(*a__[],n)
Sort the first n elements of an array. The sorting algorithm is a fast, non-recursive, run-finding mergesort.
asort0(*a__[])
Sort all elements of an array up to, but not including, the first zero (zero-terminated array).
asortb(*a__[],x)
Sort all elements of an array up to, but not including, the first instance of x.
asortl(*a__[])
Sort all elements of an array whose length is stored as its first element.
asortr(*a__[],start,end)
Sort the given range of an array.
asortr_old(*a__[],start,end)
An older version of the above which implements a recursive, non-run-finding mergesort. It is slower than the main asort...() functions but was used as a benchmark and known-correct test platform during the creation of the faster algorithm.
auniq(*v__[], a__[],n)
Fills the first array with the values from the first n values of the second array, removing any adjacent duplicate values. Returns the number of items in the first array. Note that this will be less than or equal to n. Works best on a sorted array
auniq0(*v__[], a__[])
Fills the first array with the values from the second array, which is zero-terminated, removing any adjacent duplicate values. The first array is then terminated with a zero. Returns the number of items in the first array. Works best on a sorted array
auniqb(*v__[], a__[],x)
Fills the first array with the values from the second array, which is terminated by the value given in x, removing any adjacent duplicate values. The first array is then terminated with the same terminator value. Returns the number of items in the first array. Works best on a sorted array
auniql(*v__[], a__[])
Fills the first array with the values from the second array, whose length is given in its first element. The first array is then set to be in the same format as the second, having its length as its first element. Returns the number of elements in the first array, which is the same as the aforementioned length. Works best on a sorted array
auniqr(*v__[], a__[],start,end)
Fills the first array with the values from the values found in the second array between the start and end indexes inclusive. Returns the number of elements in the first array. Works best on a sorted array
cf.bc
Rewritten for June 2012
A suite of functions for basic continued fraction analysis. Uses the undocumented pass-by-reference feature of arrays where necessary.
cf_new_(near)
Internal function; Main engine for creating continued fractions. Used by the below.
cf_tidy_()
Internal function; Main finalisation routine for tidying up after cf_new_() and others
cf_new(*cf__[],x)
Puts the continued fraction of x into the given array. For positive x, terms are always positive.
cfn_new(*cf__[],x)
Puts a better approximation continued fraction of x into the given array. Terms are positive or negative, depending on which of these obtains a closer approximation to x at each step. In this library, this form of continued fraction is referred to as a CFN (Continued Fraction Nearer).
cf_copy(*cfnew__[],cf__[])
Make a copy of the right hand array in the left, with respect to the storage rules used by this library. Some metadata may be stored in the CF arrays which ordinary array copies may lose.
cf_toggle(*cf__[],cf2__[])
Converts one type of CF into the other, i.e. CF into CFN and vice-versa, or in other words, converts the type of array created by cf_new() into the type of array created by cfn_new(). Warning: Some accuracy may be lost.
cf_value(cf__[])
Convert the continued fraction in the given array into a number
frac_from_cf(*f__[],cf__[],improper)
Fills the array f with the denominator and numerator (in that unusual order) as element 0 and element 1. If the improper flag is zero then element 2 of the array is set to be the integer part. e.g. the array would contain {2,1,5} if the specified continued fraction contained a representation of 51/2.
cf_get_convergent(cf__[],c)
Returns the c-th convergent of the given continued fraction. e.g. if the continued fraction array contains a representation of an irrational number, these convergents will represent particular approximations to that number.
get_convergent(x,c)
Finds the c-th convergent of the number x. This bypasses the need to generate a continued fraction array, but therefore can have the consequence of having to recalculate the continued fraction of x every time this function is called.
upscale_rational(x)
If x contains a rational value that was calculated with a less accurate scale, this function attempts to bring that value to full accuracy in the current scale.
cf_print(cf__[])
Displays the contents of the given array in continued fraction style. When the array contains what appears to be a representation of a quadratic irrational, will insert semicolons rather than commas after the pattern repeats, or before any pattern begins for the first time. If the global variable cf_shortsurd_ is set to non-zero this function will truncate the output of the aforementioned quadratic irrationals at the end of the first occurrence of the otherwise repeating pattern.
NB: Quadratic irrational detection is not done directly by this function, and metadata will need to be added to third-party arrays for this to work. This library adds its own metadata when putting continued fractions into arrays.
Will display an ellipsis at the end of the output if the CF is believed to be of infinite length (or at least of length beyond the current accuracy as set by scale).
cf_print_convergent(cf__[],c)
Prints a fraction representing the c-th convergent of the given array. If c is negative, will print all convergents upto and including the c-th.
cf_print_frac(cf__[],improper)
Prints a fraction representing the value of the continued fraction stored in the given array. When the improper flag is non-zero, the value displayed will, where necessary, take the form of a 'top-heavy' or improper fraction.
print_as_cf(x)
Print x as a continued fraction. This saves on creating an intermediate array if only this one output is required. Consider using other functions if other work on the continued fraction of x is required.
print_as_cfn(x)
As above, but using the CFN 'nearest sign' continued fraction format.
print_convergent(x,c)
Prints the c-th approximation to the value of x as a rational number or fraction If c is negative, prints all approximations up to the c-th.
print_frac(x,improper)
Displays x as a fraction (rational number). As in other functions the improper flag specifies whether the number is to be an improper fraction or not.
print_frac2(x,improper)
As above, but interally uses a CFN. This has the interesting side-effect of causing proper fractions with integer parts and a fractional part greater than 1/2 to be displayed with positive integer part and negative fractional part. e.g. 1 2/3 would display as 2-1/3.
cf_engel.bc
New for June 2012
Functions for generating, printing and interpreting Engel Expansions in a similar manner to ordinary continued fractions. Sister library to cf_sylvester.bc
engel_new_(mode)
Internal function used by the following three functions.
engel_new(*en__[],x)
Generate a classic EE in the given array with terms the same sign as x.
engelfall_new(*en__[],x)
Generate a secondary EE in the given array with all terms except the first having the opposite sign to x. This results in the implied Egyptian fraction having alternating signs in the sum.
engelalt_new(*en__[],x)
Generate a third kind of EE in the given array where terms alternate in sign within the EE itself. This results in the implied Egyptian fraction having pairwise alternating signs in the sum.
engel_value(en__[])
Interpret the given array as an Engel expansion and return the value that it represents.
engel_print(en__[])
Displays the contents of the given array in Engel expansion style.
Will display an ellipsis at the end of the output if the CF is believed to be of infinite length (or at least of length beyond the current accuracy as set by scale).
print_as_engel(x)
Displays x's form as a classic or primary Engel expansion.
print_as_engelfall(x)
Displays x's form as a secondary Engel expansion.
print_as_engelalt(x)
Displays x's form as the third kind of EE described above.
cf_misc.bc
New for June 2012
Miscellaneous functions that would otherwise be in the main cf.bc library.
cf_value_abs(cf__[])
Returns the value of the continued fraction stored in the array as if all terms in the continued fraction were positive.
cf_value_abs1(cf__[])
As above, but reduces all terms by 1 after making them positive. This function is written with CFNs in mind, which never have a term whose absolute value is less than 2. Will complain if a term is found to have become zero.
cf_abs_terms(*cf__[])
Irrevocably sets all terms in the given continued fraction array to their absolute value.
cf_abs1_terms(*cf__[])
Irrevocably sets all terms in the given continued fraction array to their absolute value, less one. Will complain if a term is found to have become zero.
cfn_flip_abs(x)
Transforms x through the CFN and absolute algorithms to return an alternate value.
cfn_flip_abs1(x)
As above but subtracts 1 from the absolute terms. Since this uses the CFN algorithm, no error should result.
cf_conway_box(*cf__[],x)
The inverse of the below Minkowski question mark function, this is the Conway box function □(x). This function stores the lengths of the runs of 1s and 0s into the given array and formats said array as if it were the continued fraction of another number instead.
conway_box(x)
Transforms a number using the above.
cf_question_mark(cf__[])
The Minkowski question mark function ?(x). Transforms a continued fraction into a binary string by treating the terms as lengths of runs of 1s or 0s in a new number. Since quadratic irrationals have repeating continued fraction expansions, their question-mark transformations have repeating binary expansions and are therefore necessarily transformed into rational numbers.
question_mark(x)
As above but transforms a number rather than parsing a continued fraction.
cf_sylvester.bc
New for June 2012
Functions for generating, printing and interpreting Sylvester Expansions in a similar manner to ordinary continued fractions. Sister library to cf_engel.bc
sylvester_new_(mode)
Internal function used by the two following
sylvester_new(*sy__[],x)
Create a greedy egyptian fraction of x also known as a Sylvester expansion in the given array.
sylvester2_new(*sy__[],x)
As above but creates a secondary kind of SE where all terms after the first are of opposite sign to x.
sylvester_value(sy__[])
Interpret the given array as a Sylvester expansion and return the value that it represents.
sylvester_print(sy__[])
Displays the contents of the given array in an easily readable style.
Will display an ellipsis at the end of the output if the CF is believed to be of infinite length (or at least of length beyond the current accuracy as set by scale).
print_as_sylvester(x)
Displays x's form as a classic or primary Sylvester expansion.
print_as_sylvester2(x)
Displays x's form as a secondary Sylvester expansion.
collatz.bc
Rewritten for June 2012
A suite of functions for very basic experimentation with the Collatz conjecture.
All functions here use the global variable collatz_mode_ to determine which Collatz ruleset is to be used.
By default this is set to 1 and the rules are the standard Collatz rules of:
even x → x/2, odd x → 3x+1.
When set to 2, the rules become the condensed rules:
even x → x/2, odd x → (3x+1)/2.
When set to zero, the rules become:
even x → oddpart(x), odd x → oddpart(3x+1), where oddpart is what remains of a number when it is divided by 2 until it cannot be divided further.
collatz_next(x)
Returns the next hailstone on from x
collatz_prev(x)
Returns one of the possible previous hailstones from x, choosing the lower of the available options.
is_collatz(x)
If the Collatz conjecture is true, this function will return 1 (true) for all positive integers. Will return -1 (true but negative) if a negative number reaches -1, though many negative numbers reach a loop and so, where detected, this function may return 0 (false)!
collatz_print(x)
Displays all hailstones from x down to 1 or -1, or until a loop is found. NB: This function fits with the system used elsewhere in these files and so should have its value assigned to a variable in order to avoid printing its return value, which is 0.
collatz_root(x)
Returns the lowest number reached by the Collatz iteration whether that be 1, -1 or the smallest member of the loop. Note that all positive integers are thought to reach 1.
collatz_loopsize(x)
Returns the size of the loop that that Collatz iteration eventually becomes caught within when starting from x. For positive integers under the standard Collatz iteration, this loopsize is 3 because of the 4 → 2 → 1 → 4 cycle.
collatz_chainlength(x)
Returns the number of iterations before reaching a loop or a terminating condition.
collatz_magnitude(x)
Returns the highest point reached during the Collatz iteration.
collatz_sum(x)
A hold-over from a previous version of this library. Returns the sum of all hailstones reached during the iteration.
is_collatz_sg(x)
A combination of all of the above functions, setting global variables by same names. In the case of collatz_print, this must be set by the user, zero for disabled, and non-zero if this function is to print all hailstones on the way to the terminating condition. Returns the same value as would is_collatz()
collatz_continuous.bc
New for June 2012
Most functions in this library have two versions; One which relies on the global collatz_mode_ variable as described in collatz.bc and another which has a secondary parameter for otherwise specifying the power of two of the divisor (if any) on the odd step. See notes in code for more details.
Each of the non-inverse functions is designed to provide a continuous version of the single Collatz iteration function. This necessarily means that each of mappings has fixed points at certain non-integer values of their first parameter. i.e. x = f(x)
collatz_arccos(y)
Inverse of collatz_cos(x)
collatz_arccos_(y,k)
Inverse of collatz_cos_(x,k)
collatz_arcpcos(y)
Inverse of collatz_pcos(x)
collatz_arcpcos_(y,k)
Inverse of collatz_pcos_(x,k)
collatz_cos(x)
A continuous version of the single Collatz iteration function using cosine interpolation to obtain values for non integer x.
collatz_cos_(x,k)
As above but with more control over the divisor on the odd step. k may take any value, integer or otherwise.
collatz_pcos(x)
A continuous version of the single Collatz iteration function using a clumsier, piecewise cosine interpolation to obtain values for non-integer x.
collatz_pcos_(x,k)
As above but with more control over the divisor on the odd step. k may take any value, integer or otherwise.
collatz_invlin(y)
Inverse of collatz_lin(x)
collatz_invlin_(y,k)
Inverse of collatz_lin_(x,k). Is used by all other inverse functions as a starting point.
collatz_invlinb(y)
Inverse of collatz_linb(x)
collatz_invlinb_(y,k)
Inverse of collatz_linb_(x,k)
collatz_lin(x)
Piecewise linear interpolation of the single Collatz iteration function allowing for the generation of 'iterations' for non-integer x. Unlike for cosine interpolation, piecewise is a better choice for linear interpolation.
collatz_lin_(x,k)
As above but with more control over the divisor on the odd step. k may take any value, integer or otherwise.
collatz_linb(x)
An alternative, non-piecewise, pseudo-linear interpolation of the single Collatz iteration function. While this has an arguably beautiful closed form, much like the non-piecewise cosine interpolation, this function, when graphed, does not follow a straight path between the integer values either side of x.
collatz_linb_(x,k)
As above but with more control over the divisor on the odd step. k may take any value, integer or otherwise.
collatz_piecewise__(y,k)
Internal function; Used by the inverse functions for solving piecewise transformations.
complex.bc
A second attempt at creating and working with complex numbers in bc. Uses arrays and the undocumented pass-by-reference feature to store the real and imaginary parts of a number. The downside to this method is that the syntax is somewhat unwieldy: When a complex return value is required, the first parameter is always *c__[] and the return value is stored within the first two elements of the supplied array, rather than being returned in the usual bc way.

Some 'constants' are predefined by this library; These are complex0[], complex1[], complex2[], complexi[], complexomega[] and complexomega2[]. These are zero, one, two, imaginary unit i and roots of positive unity ω and ω2.

arctan2(x,y)
Two-parameter inverse tangent; Takes two ordinary numbers and returns an ordinary number. This is used by a few of the below functions.
cabs(c__[])
Returns an ordinary number which is the absolute value or magnitude of the complex number stored in the first two elements of the given array.
cadd(*c__[],a__[],b__[])
Stores the sum of the second and third parameters into the first; i.e. c = a + b
caddassign(*c__[],b__[])
Adds the value of the second parameter to the value already stored in the first; i.e. c += b
carccis(*c__[],x__[])
Inverse of ccis(); If x is of form cos(c)+i.sin(c), calculate a possible value for c.
carccos(*c__[],x__[])
Stores the inverse cosine of x in c.
carccosec(*c__[],x__[])
Stores the inverse cosecant of x in c.
carccosech(*c__[],x__[])
Stores the inverse hyperbolic cosecant of x in c.
carccosh(*c__[],x__[])
Stores the inverse hyperbolic cosine of x in c.
carccotan(*c__[],x__[])
Stores the inverse cotangent of x in c.
carccotanh(*c__[],x__[])
Stores the inverse hyperbolic cotangent of x in c.
carcsec(*c__[],x__[])
Stores the inverse secant of x in c.
carcsech(*c__[],x__[])
Stores the inverse hyperbolic secant of x in c.
carcsin(*c__[],x__[])
Stores the inverse sine of x in c.
carcsinh(*c__[],x__[])
Stores the inverse hyperbolic sine of x in c.
carctan(*c__[],x__[])
Stores the inverse tangent of x in c.
carctanh(*c__[],x__[])
Stores the inverse hyperbolic tangent of x in c.
carg(c__[])
Returns an ordinary number which is the angle or argument of the complex number in c.
carrayget(*c__[],a__[],i)
Rudimentary complex array handling; Second parameter is treated as an array of complex numbers, which in bc terms is an array of numbers that for each pair of elements, a complex number is represented. This function stores the i-th complex entry in the array into c. i.e. c = a[i]
carrayset(*a__[],i,c__[])
Rudimentary complex array handling; First parameter is treated as an array of complex numbers, which in bc terms is an array of numbers that for each pair of elements, a complex number is represented. This function stores the complex number found in c into the i-th complex entry of that array. i.e. a[i] = c
cassign(*c__[],x__[])
Stores a copy of the first two elements of x into the first two elements of c. i.e. set one complex number equal to another.
ccis(*c__[],x__[])
Stores cos(x)+i.sin(x) into c. To calculate the same for x as an ordinary bc number, see the cis() function.
cconj(*c__[],x__[])
Stores the complex conjugate of x into c.
cconjself(*c__[])
Turn a complex number into its conjugate.
ccos(*c__[],x__[])
Stores the cosine of x in c.
ccosec(*c__[],x__[])
Stores the cosecant of x in c.
ccosech(*c__[],x__[])
Stores the hyperbolic cosecant of x in c.
ccosh(*c__[],x__[])
Stores the hyperbolic cosine of x in c.
ccotan(*c__[],x__[])
Stores the cotangent of x in c.
ccotanh(*c__[],x__[])
Stores the hyperbolic cotangent of x in c.
cdiv(*c__[],a__[],b__[])
Stores the result of dividing the second parameter by the third into the first; i.e. c = a / b
cdivassign(*c__[],b__[])
Divides the value stored in the first parameter by the value in the second; i.e. c /= b
cequal(a__[],b__[])
Boolean function; Returns 1 if the complex numbers represented by the parameters are equal.
cexp(*c__[],x__[])
Stores the x-th power of e into c. Exponential.
cexpself(*c__[])
Raises e to the power of the parameter and then stores it back there.
cintpow(*c__[], z[], n)
Stores the nth power of z in c; i.e. c = z^n (where n is integer). This is faster than using the cpow() function.
cintpowassign(*c__[],n)
Raises c to the nth integer power; i.e. c = c^n This is faster than using the cpowassign() function.
cis(*c__[],x)
Stores cos(x)+i.sin(x) into c. To calculate the same for x as a complex number, see the ccis() function.
cln(*c__[],x__[])
Stores the natural logarithm of x in c.
clnself(*c__[])
Takes the natural logarithm of c and replaces the original value with it.
cmul(*c__[],a__[],b__[])
Stores the product of the second and third parameters into the first; i.e. c = a.b
cmulassign(*c__[],b__[])
Multiplies the value stored in the first parameter by the value in the second; i.e. c = c.b
cneg(*c__[],x__[])
Stores the value of -1 times x in c.
cnegself(*c__[])
Inverts the sign of the given complex number.
cpow(*c__[],a__[],b__[])
Stores the result of raising the second parameter to the power of the third into the first; i.e. c = a^b
cpowassign(*c__[],b__[])
Raises the first parameter to the power of the second; i.e. c = c^b
csec(*c__[],x__[])
Stores the secant of x in c.
csech(*c__[],x__[])
Stores the hyperbolic secant of x in c.
csgn(*c__[],x__[])
Stores the complex sign, or unit circle intersection point of x in c.
csgnself(*x__[])
As above but performs the operation on the given number and overwrites.
csin(*c__[],x__[])
Stores the sine of x in c.
csinh(*c__[],x__[])
Stores the hyperbolic sine of x in c.
csqrt(*c__[],x__[])
Stores the principle square root of x in c.
csqrtself(*c__[])
Overwrites the parameter with its principle square root.
csquare(*c__[],x__[])
Stores the square of x in c. Faster than the ...pow...() functions or multiplying.
csquareself(*c__[])
Overwrites the parameter with its square.
csub(*c__[],a__[],b__[])
Stores the result of subtracting the third parameter from the second into the first; i.e. c = a - b
csubassign(*c__[],b__[])
Subtracts the values of the second parameter from the value held in the first; i.e. c = c - b
ctan(*c__[],x__[])
Stores the tangent of x in c.
ctanh(*c__[],x__[])
Stores the hyperbolic tangent of x in c.
imag(c__[])
Returns, as a standard bc number, the imaginary part of c.
int(n)
Returns the integer part of an ordinary bc number.
makecomplex(*c__[],r,i)
Creates a complex number in c from real and imaginary parameters r and i, repectively.
makeomega()
Technically an internal function. Sets the values of intentionally predefined array-based complex numbers. These are complexomega[] and complexomega2[] which are the complex third roots of unity. This function should be called if these two arrays are corrupted for some reason.
mod(n,m)
Returns n modulo m. Return value and parameters are all standard bc numbers, and indeed should all be integers.
printc(c__[])
Prints the complex number stored in c in a human readable format.
real(c__[])
Returns, as a standard bc number, the real part of c.
cosconst.bc
The cosine constant to a large number of decimal places, where x = cos(x)
None
There is only a constant called cosconst
digits.bc
Minor update for February 2013
Treat numbers as strings of digits. Some of the definitions below are not in strict alphabetical order. This is so that concepts are introduced in a more logical order.

Many functions will operate with bijective (zero-less) number bases, which can be activated through use of the new global bijective flag. By default this is zero. Set to non-zero to alter the behaviour of those functions which support it.

Some functions operate correctly with negative number bases, e.g. negabinary.

base_check_()
Internal function; used by most functions to check the sanity of the base parameter which most of these have.
bmod_(x,y)
Internal function; Returns x modulo y and sets global variable bdiv_ to the value of the division which found the modulus. In bijective mode the return range is 1..y rather than 0..y-1, and bdiv_ is one less than the true division when the modulus is equal to y. For positive x and y, x == bmod_(x,y)+y*bdiv_ is always true.
cantor(basefrom, baseto, x)
Treat x's representation in basefrom as a representation in baseto and return the resulting number, i.e. reinterpret the number.
Will always convert successfully to a larger base, but the reverse is often not possible, and a warning will result when data loss occurs.
Warnings can be turned off by setting the global cantorwarn_ variable to 0. It is set to 1 by default.
Some warnings can be mitigated through use of the similar global cantormod_ variable. If this is set to 1, calculated digits are truncated, modulo baseto, meaning some information will be lost, but no left-carries will occur.
x can be interpreted as if basefrom is a bijective base (see output_formatting.bc for ways to display numbers as bijective) if global variable bijective is set to 1. Note that a warning will result if x is non-integer in bijective mode; there is no correct way to interpret such a value.
cantor_trans(basefrom, baseto, mul, cons, x)
As above, but an extra linear transform can be applied to the digits as they are converted between the bases. They will be multiplied by the given multiplier and then added to the given constant.
cantor_trans_(d)
Internal function; Acts as an error checker / reporter for the above.
cantor_trans_negabase_(basefrom, baseto, mul, cons, x)
Internal function; Used by cantor_trans to be able to handle negative integer bases. Note that the parameters match.
digit_sum(base,x)
Add the digits of x when interpreted in the specified base. Repeated applications of this function would derive the "digital root".
digit_product(base,x)
As above, but multiply rather than add. e.g. 235 -> 2*3*5 = 30 in base ten. In non-bijective bases this function returns zero much of the time since there is a strong chance that one or more of the digits in the number is a zero. See digits_misc.bc for two alternatives to this function which do not have this disadvantage.
digit_distance_(base,x,y,sh)
Internal function: Engine used by both of the below
digit_distance(base,x,y)
Adds the differences of respective digits of x and y in the given base and returns the result. e.g. 246(-)176 ->|2-1|+|4-7|+|6-6| = 1+3+0 = 4 in base ten. This can be considered an extension to the concept of Hamming distance, which merely counts the number of differences. See also base_hamming() in logic_otherbase.bc and hamming() in logic.bc for pure difference functions.
digit_sdistance(base,x,y)
In modular arithmetic, there it can be argued that there are often two possible positive solutions to |x-y|; One of these matches the usual definition, but the other has a value of modulus-|x-y|. e.g. modulo ten, |1-8|=7, but wrapping around zero, we can move 8 to 9, 9 to 0 and 0 to 1, which is only three steps.
This function therefore, is the same as the above but uses the smaller of the two options in the summation for the distance between digits.
digits(base,x)
Find the number of digits in x's representation in the given base.
int_catenate(base, x,y)
Splice two integer representations together in the specified base so that x is before y.
int_left(base, x, count)
Returns the leftmost digits of x in the given base, specified by the given count.
int_mid(base, x, start, end)
Returns digits of x in the given base, counting in from the left, starting and ending at the given digit positions.
int_right(base, x, count)
Returns the rightmost digits of x in the given base, specified by the given count.
is_palindrome(base,x)
Determine whether x reads the same forwards and backwards in the given base
is_pseudopalindrome(base,x)
Determine whether x reads the same forwards and backwards in the given base, or could read the same each way if zeroes were prepended to the number (which wouldn't actually change its value).
is_substring(base,large,small)
Determine whether the digits of the smaller number appear, in order, within the digits of the larger number, all in the given base.
make_even_palindrome(base, x)
Turn x into a unique palindrome with an even number of digits in the given base.
make_odd_palindrome(base, x)
Turn x into a unique palindrome with an odd number of digits in the given base.
map_palindrome(base, x)
Generate a unique palindrome from x in the given base. This function maps the integers onto the palindromes on a one-to-one basis.
reverse(base,x)
Reverse the digits of x in the current base. Zeroes at the end of x will be lost.
unmap_palindrome(base, x)
Inverse function of map_palindrome; Maps the domain of palindromes in the given base back into the integers.
digits_describe.bc
describe_(opt,base,x)
Internal function for both modes of describing numbers
describe_countfirst(base,x)
Generates a number (in the specified base) which describes x by putting the digit count of each digit of x before the actual digit of x. This is the standard, well known look-and-say sequence. e.g. 111 -> 31 (three 1s); 1123 -> 211213 (two 1s, one 2, one 3). A warning will result if a digit count is too large for the specified base.
describe_countlast(base,x)
Generates a number (in the specified base) which describes x by putting the digit count of each digit of x after the actual digit of x. This is the alternative look-and-say sequence. e.g. 111 -> 13 (1, three times); 1123 -> 122131 (1 twice, 2 once, 3 once) Again, a warning will result if a digit count is too large for the specified base.
parserle_(opt,base,x)
Internal function for both modes of interpreting the above description numbers. The name comes from "Parse RLE" or Parse Run Length Encoding. The irony of the function name being hard to read (parse) has been left uncorrected as it is amusing to the author as well as the same length (in letters) as "describe".
parserle_countfirst(base,x)
Inverse of describe_countfirst(); Interprets the value in x as a description (in the specified base) of a number, which is calculated and returned. A warning will result if x is not interpretable.
parserle_countlast(base,x)
Inverse of describe_countlast(base,x); Interprets the value in x as a description (in the specified base) of a number, which is calculated and returned. A warning will result if x is not interpretable.
digits_happy.bc
Very minor update for February 2012
A suite of functions for investigating the so-called Happy Numbers. This library supports the bijective global variable used elsewhere in these files. When this variable is non-zero (and thus bijective mode is enabled), numbers are also considered to be happy if the iteration reaches the base and not just 1.
base_check_happy_()
Internal function; used by most functions to check the sanity of the base parameter which most of these have.
is_happy(base,pow,x)
Returns 1 (true) if x is happy in the given base when each digit is raised to the given power, 0 (false) otherwise. The original definition of happiness involves base ten and a power of two (squaring).
happy_chainlength(base,pow,x)
Returns the number of iterations required to reach happiness for the given x. Will return a negative number of iterations if x is unhappy, which is the number of iterations required before the algorithm "realises" that the number is unhappy.
happy_loopsize(base,pow,x)
Should probably be called unhappy_loopsize, since this will return the number of iterations in the loop encountered by numbers that are not happy.
happy_print(base,pow,x)
Shows all the iterations down to 1 if the number is happy, or else stops once a loop is detected.
happy_root(base,pow,x)
Shows 1 if the number is happy, or else the smallest number in the loop that an unhappy number becomes trapped within.
is_happy_sg(base,pow,x)
sg = "set globals": This function is all of the above rolled into one, and will set global variables by the names of the above functions (e.g. happy_root) for the parameters given. If the global variable happy_print is set to 1, then this function will also behave as the happy_print() function and display the value of the iterations. Set to 0 to turn the feature off again. Like is_happy(), returns 1 if x is happy and 0 otherwise.
digits_misc.bc
Very minor update for February 2012
Some of the more obscure digit stringification functions you may wish to encounter, removed from digits.bc as interesting but unnecessary, or rolled in from various old bc files long removed from this page.
base_check_misc_()
Internal function; used by most functions to check the sanity of the base parameter which most of these have.
append_all(base,x)
The digit string equivalent of the triangular numbers or the factorials. Appends all representations of the numbers from 1 to x in the current base to each other. e.g. assuming base ten, append_all(10, 15) = 123456789101112131415
calcsegments(base,x)
Returns the number of segments that would be 'lit' on a seven-segment-per-number calculator display. Customised to support bases up as far as 36, although no calculator goes any further than 16. Adds one for the negative sign since all calculators need a segment to show that.
count_digit(base,x,digit)
Returns the number of occurrences of a given digit within an integer, x, in the given base. e.g. 52726620 has 3 occurrences of the digit 2 in base ten.
count_digits(*d__[],base,x)
Uses the undocumented pass-by-reference for arrays to store a count of all types of digit in the given base that can be found in an integer, x. e.g. 10110101010010101 in binary would result in d containing {8,9} because there are eight zeroes and nine ones.
digit_product1(base,x)
An alternative to digit_product() in digits.bc. Rather than multiply the digits immediately, one is added to each before the multiplication and then one is subtracted from the final product, ensuring a non-zero result. e.g. 235 -> (2+1)(3+1)(5+1)-1 = 3*4*6 - 1 = 71 in base ten.
digit_product2(base,x)
Another alternative to digit_product() and the above. All digits are translated into their corresponding odd number, multiplied and then mapped back from the odd integers to the integers again. e.g. 13462 -> ( (2*1+1)(2*3+1)(2*4+1)(2*6+1)(2*2+1)-1 )/2 = (3*7*9*13*5 - 1)/2 = 6142 in base ten.
is_negapalindrome(base,x)
Determine whether the opposing pairs of digits, (counted in from either end) sum to one less than the given base. e.g. 147258 is a negapalindrome in base ten since 1+8 = 4+5 = 7+2 = 9 = 10 - 1
is_pseudonegapalindrome(base,x)
Determine whether x is a negapalindrome in the given base should any number of zeroes are prepended to the number. These would tie in with any digits one less than the base found at the end of x, and wouldn't change x's value.
is_negapalindrome2(base,x)
Alternate definition of negapalindrome, where opposing pairs of digits must sum to the base itself, rather than one less.
map_negapalindrome(base, x)
Generate a unique negapalindrome from x in the given base. This function maps the integers onto the negapalindromes on a one-to-one basis.
pdhi(x)
Pan digital halving index. Determine the number of times that x must be halved before its decimal expansion (or the expansion in the base specified by ibase) contains all possible digits. Warning: May hang for some values of x
pdmi(x,m)
Pan digital multiplying index: Determine how many times x must be multiplied by m before it is pandigital. Warning: May hang for some values of x and m
sdp(base,x)
Swap Digit Pairs: Takes the representation of x in the given base and reverses every pair of digits.
sort_digits_asc(base,x)
Sort the digits of x into ascending order in the given base. Zeroes at the end of x will be lost.
sort_digits_desc(base,x)
Sort the digits of x into descending order in the given base.
split_digits(*d__[],base,x)
Using the given base, store x into the given array in an unambiguous manner which does not lose any information. For more details, see the source code.
join_digits(d__[])
Returns the number which has been stored into an array by the split_digits function. Basimal point and original base are store within the array hence only one parameter being required.
stripbm1s_(base,x)
Internal function used by the negapalindrome family.
unmap_negapalindrome(base, x)
Inverse function of map_negapalindrome; Maps the domain of negapalindromes in the given base back into the integers.
factorial.bc
Updated for February 2013
A suite of functions for calculations involving the factorial and its relatives. Migrated from funcs.bc and then expanded upon. Still requires funcs.bc to work correctly.
factorial(x)
An approximation to the factorial function over the reals. Is accurate as possible for all integers and half-integers, but interpolates otherwise. Accuracy versus speed can be balanced by changing the value of the global factorial_substrate_ variable. Set to 0, will use the cheapest approximation for interpolation. Values of 1, 2 (default) and 3, are increasingly slower but more accurate. A value of 4 (the maximum) will ensure calculations take as long as necessary to find an accurate answer.
lncombination(n,r)
Calculates the logarithm of the combination() function using the lnfactorial() function so that larger values of n and r may be calculated for more quickly.
lnpermutation(n,r)
Calculates the logarithm of the permutation() function using the lnfactorial() function so that larger values of n and r may be calculated for more quickly.
lnfactorial(x)
Calculates the logarithm of the factorial() function in a way generally faster than ln(factorial(x)), but with the same caveats as before: Is accurate as possible for all integers and half-integers, but interpolates otherwise.
catalan(n)
Returns the nth Catalan number
combination(n,r), int_combination(n,r)
Calculates the binomial coefficient nCr. i.e. How many ways can r objects be chosen from n objects without regard to order? The non-integer function is slower but uses the factorial function to a closely approximated calculation for non integral parameters.
double_factorial(x)
Calculates the so-called double factorial (x!! = x.(x-2).(x-4)..{2 or 1}) with the same caveats as for the factorial() and other functions: Is accurate as possible for all integers and half-integers, but interpolates otherwise.
eulergamma()
Gives Euler's γ (0.5772156...) to the number of decimal places specified by the current scale. Uses global variable eulergamma_ to cache the value to save on recalculating every time this is called. The variable always contains a value of accuracy greater or equal to the current scale. Caution: this function will be very slow on first running as it calculates the value in real time.
factorial_substrate_(n), lnfactorial_substrate_(n)
Internal functions. These derive an appropriate approximation specified by the aforementioned global factorial_substrate_ variable, and return the result.
fast_inverse_factorial(x)
A very approximate inverse to the factorial() function. Developed from an idea by David W. Cantrell.
fast_inverse_lnfactorial(x)
A very approximate inverse for the lnfactorial() function.
gfactorial(n)
A rough, quick and dirty approximate to the factorial function using the below. [Substrate 0]
gosper(x)
Gosper's approximation to the natural logarithm of the factorial function. [Substrate 0]
int_multifactorial(y,x)
Quick and dirty function to calculate the y'th multifactorial of x.
inverse_factorial(f)
Uses and improves on the result given by fast_inverse_factorial(x) to give a an almost exact inverse. Much slower than the latter, but useful at high substrates for providing as accurate an answer as possible.
inverse_lnfactorial(x)
Uses and improves on the result given by fast_inverse_lnfactorial(x) to give a an almost exact inverse. Much slower than the latter, but useful at high substrates for providing as accurate an answer as possible.
lcmultorial(x)
Calculate the lowest common multiple of all integers from 1 to x. This is an integer-only function. For an attempt at a more continuous function, see lcmultorialc() in orialc.bc.
nemes(x)
Gergo Nemes' excellent approximation to the natural logarithm of the factorial function. [Substrate 1]
nemfactorial(n)
Uses the above to calculate an approximation to the factorial function. [Substrate 1]
permutation(n,r), int_permutation(n,r)
How many ways can r objects be chosen from n objects when the order of choosing is important? The non-integer function is slower but uses the factorial function to a closely approximated calculation for non integral parameters.
spouge(n)
John L. Spouge's approximation to the natural logarithm of the Gamma function, adjusted to be an lnfactorial approximation; This version calculates to within scale significant figures. [Substrate 3]
spougefactorial(n)
John L. Spouge's approximation to the Gamma function, adjusted to be a factorial approximation; This version calculates to within scale significant figures. [Substrate 3]
spougefactorialx(n)
John L. Spouge's approximation to the Gamma function, adjusted to be a factorial approximation; This version calculates to within scale decimal places, which may take much longer than the x-less counterpart above. [Substrate 4]
spougex(n)
John L. Spouge's approximation to the natural logarithm of the Gamma function, adjusted to be an lnfactorial approximation; This version calculates to within scale decimal places, which may take much longer than the x-less counterpart above. [Substrate 4]
spouge_(n,l,exact)
Internal function used by the other Spouge functions. [Substrate >3]
stielfactorial(n)
Thomas J. Stieltjes' approximation to the Gamma function, adjusted to be a factorial approximation. [Substrate 2]
stieltjes(n)
Thomas J. Stieltjes' approximation to the natural logarithm of the Gamma function, adjusted to be an lnfactorial approximation. [Substrate 2]
subfactorial(n)
Subfactorial / Derangement counting function; Gives how many ways a number of items can be rearranged such that no item is in its original place
factorial_gamma.bc
New for February 2013
Near-aliases for some functions in factorial.bc, giving the functions more classical function names.
gamma(x)
An approximation to the Gamma function over the reals. Is accurate as possible for all integers and half-integers, but interpolates otherwise. Accuracy versus speed can be balanced by changing the value of the global
lngamma(x)
Calculates the logarithm of the gamma() function in a way generally faster than ln(gamma(x)), but with the same caveats as before: Is accurate as possible for all integers and half-integers, but interpolates otherwise.
inverse_gamma(f)
Relatively slow but accurate calculation of the inverse of the Gamma function for positive real values.
inverse_lngamma(x)
Relatively slow but accurate calculation of the inverse of the lnGamma function for positive real values.
beta(x,y)
Euler's Beta function
lnbeta(x,y)
Calculates the logarithm of the beta() function in a way generally faster than ln(beta(x,y))
fibonacci.bc
Fibonacci and Lucas functions; Originally part of funcs.bc Still requires funcs.bc to work correctly.
fibonacci(n)
An extension of the fibonacci numbers over the reals without stepping into complex numbers, which would be necessary for a true extension.
inverse_fibonacci(f)
An inverse to the fibonacci function. Provides incorrect answers for non integer values of f when f < 1.
inverse_lucas(l)
An inverse to the lucas function. Provides incorrect answers for non integer values of l when l < 2.
lucas(n)
Returns the n'th Lucas number. Continuous over the reals, like its cousin the fibonacci function
funcs.bc
Minor update for February 2013
A large suite of functions to complement the bc standard library. Unlike the standard library (activated with bc -l), all function names are spelled out in full. Full name aliases for the standard library functions are provided.
abs(x)
Absolute value of a number
arccos(x)
Inverse cosine
arccosec(x)
Inverse cosecant
arccosech(x)
Inverse hyperbolic cosecant
arccosh(x)
Inverse hyperbolic cosine
arccotan(x)
Inverse cotangent (single variable)
arccotan2(x,y)
Inverse cotangent (two axes)
arccoth(x)
Inverse hyperbolic cotangent
arcgudermann(x)
Inverse of the Gudermann function
arcsec(x)
Inverse secant
arcsech(x)
Inverse hyperbolic secant
arcsin(x)
Inverse sine
arcsinh(x)
Inverse hyperbolic sine
arctan(x)
Inverse tangent (single variable). This is an alias for bc's own a() function.
arctan2(x,y)
Inverse tangent (two axes)
arctanh(x)
Inverse hyperbolic tangent
arigeomean(a,b)
Arithmetic-geometric mean
besselj(n,x)
Bessel J function. This is an alias for bc's own j() function.
ceil(x)
Ceiling function: returns the next integer greater than or equal to x
converse_poly(x,r)
converse of poly; solves poly(s,x)=r for s. i.e. if the xth polygonal number is r, how many sides has the polygon? e.g. if the 5th polygonal number is 15, converse_poly(5,15) = 3 so the polygon must have 3 sides! (15 is the 5th triangular number)
cbrt(x)
Cube root; Specific and streamlined version of root(x,3), included by popular demand. Will always find an integer cube root without loss of accuracy no matter how large is x.
cos(x), c(x)
Cosine; Uses the pi() function below in order to provide faster calculation than bc's library c() function. Overrides that same function due to the speed increase.
cosec(x)
Cosecant
cosech(x)
Hyperbolic cosecant
cosh(x)
Hyperbolic cosine
cotan(x)
Cotangent
coth(x)
Hyperbolic cotangent
evenpart(x)
Find the largest power of two within x, i.e. the even part of x
exp(x)
Exponential function ex. This is an alias for bc's own e() function.
fastfracpow_(x,y)
Internal function for faster / alternative methods of raising a number to powers between 0 and 1.
fastintpow_(x,y), fastintpow__(x,y)
Internal functions for faster / alternative methods of raising bc numbers to integer powers. See source code for more information. Used by the pow() function below.
fastintroot_(x,y)
Internal function for finding integer roots in a faster and less scale-reliant way than other internal functions. As above, is used by pow().
floor(x)
Floor function. Finds the integer less than or equal to x.
frac(x)
Finds the fractional part of number, discarding the integer part. Always returns a non-negative answer.
gcd(x,y), int_gcd(x,y)
Calculate the GCD (Greatest Common Divisor) of x and y.
gudermann(x)
The Gudermann function which links hyperbolic and common trigonometric functions.
id_frac2_(y)
Internal function. Helps determine whether the fractional part of a number is most likely odd/even, odd/odd, even/odd or irrational. See comments in code for more.
int(x)
Finds the integer part of x, always rounding towards zero. See ceil and floor for more useful functions.
inv_arigeomean(n, y)
Inverse of the arithmetic-geometric mean; Given the arithmetic-geometric mean of two numbers, n, and one of the two numbers, y, finds the value of the unknown number.
inverse_poly(s, r)
"Polygonal root": If a polygonal number with s sides has area r, how many elements are along each side? For s = 4 this is the same as the square root, and for s = 3, this is the same as the trirt function.
lambertw0(x)
The zero branch of the Lambert W function, i.e. the inverse of xex.
lambertw0_exp(x)
Faster, more advisable alternative for calculating lambertw0(exp(x)).
lambertw_1(x)
The minus one branch of the Lambert W function
lcm(x,y), int_lcm(x,y)
Calculate the LCM (Lowest/least Common Multiple) of x and y.
ln(x)
A speed improvement to bc's own l() Natural Logarithm function, though the latter is still used by this function. Complains when given unexpected values.
log(base,x), int_log(base,x)
Find the logarithm of x to the given base.
oddpart(x)
Finds the odd part of a number, removing any powers of two, e.g. the odd part of 60 is 15
phi()
Gives the golden ratio φ (1.618033...) to the number of decimal places specified by the current scale.
pi()
Gives π (3.141592...) to the number of decimal places specified by the current scale. Uses global variable pi_ to cache the value to save on recalculating every time this is called. The variable always contains a value of accuracy greater or equal to the current scale.
poly(s, x)
Return the x'th s-sided polygonal number, e.g. the 10th triangular number = poly(3,10)
pow(x,y)
Returns an extremely close approximation (completely accurate in the case of integer parameters) to xy; Copes very well with negative numbers, fractional exponents etc. always returning a real or exact integer root where possible. Will complain and return zero otherwise.
powroot(x)
Solves x = yy for y.
psi()
Gives the alternate golden ratio ψ (-0.618033...) to the number of decimal places specified by the current scale.
pyth(x,y)
Pythagoras: Calculates the hypotense of a right angled triangle whose other sides are x and y.
pyth3(x,y,z)
Pythagoras 3D: Calculates the long diagonal of a cuboid whose sides are x, y and z.
remainder(x,y), int_remainder(x,y)
Calculates the remainder when x is divided by y. The non-integer version works in a more intuitive manner than bc's built in % (modulus) operator.
root(x,y)
Returns an extremely close approximation to y√x; Copes very well with negative numbers, fractional exponents etc. always returning a real or integer root where possible. Will complain and return zero if there is a problem.
round( x,y)
Round x to the nearest multiple of y.
round_down(x,y)
Round x to the multiple of y less than or equal to x.
round_up( x,y)
Round x to the multiple of y greater than or equal to x.
sec(x)
Secant
sech(x)
Hyperbolic Secant
sgn(x)
Returns the sign of x; -1 for negative, 0 for zero, 1 for positive
sin(x)
Sine; Is an alias for bc's own s() function.
sinh(x)
Hyperbolic sine
tan(x)
Tangent
tanh(x)
Hyperbolic tangent
tet(n)
The n'th tetrahedral number
tetrt(t)
"Tetrahedral root": Given a tetrahedral number, returns its index in the sequence of tetrahedral numbers. Akin to a cube root.
tri(x)
The x'th triangular number
tri_pred(t)
Given a triangular number t, returns the next triangular number. Works also for non triangular numbers, providing a continuum.
tri_step_(t,s)
Internal function: Used by the preceding and succeeding entries here...
tri_succ(t)
Given a triangular number t, returns the previous triangular number. Works also for non triangular numbers, providing a continuum.
trirt(x)
"Triangular root": Given a triangular number, returns its index in the sequence of triangular numbers. Akin to a square root.
w(x)
In the manner of bc's own single-letter functions s(), c(), a(), l(), e() and j(), this provides access to the lambertw... functions, choosing the most logical branch; Minus one for negative x, Zero for positive and zero x.
w_e(x)
A shorthand alias for the lambertw0_exp() function, i.e. a better alternative to w(e(x)).
intdiff.bc
Perform numerical integration and differentiation of a single variable function.
f(x)
All ?fxdx functions here automatically look for a function called f to perform their operations upon. Since bc allows re-definition of functions, redefining f(x) to be an alias of the function to be used is recommended before using the other functions. e.g. define f(x){return sqrt(x)}; ifxdx(2,3)
dfxdx(x)
Return the value of the first derivative of f at x.
glai(p,q,r)
Guess Limit At Infinity: given three convergents to a limit, this function attempts to extrapolate the limit at infinity. e.g. glai(63.9, 63.99, 63.999) returns 64. Uses global variable glaitalk to comment on and warn about interesting situations. Set this to 0 to turn it off.
ifxdx(a,b)
Return the indefinite integral (i.e the area under the curve) of f between a and b. A global variable called depth is used here (akin to bc's own scale variable), which determines how deep the calculation should go. It is set at an acceptable (for 2010) value already. The user changes it at their own risk as calculation time grows exponentially in proportion to it.
ifxdx_g(a,b)
As above but uses glai to save on calculations.
interest.bc
Early in 2011, Randy Rysavy, a visitor to this site, contacted me enquiring about the possibility of adding financial functions to GNU bc. He included some of his own sample code, which gave me the impetus to produce this suite. Randy's own functions aren't included here, although there are equivalents and many more besides.
As in other places on this site, the list of functions is not entirely in alphabetical order so to introduce concepts in a more logical order.
Functions here rely fairly heavily on exponential related functions found in funcs.bc.
fraction_to_rate(fraction)
This library expects interest rates to be given to functions in the form of 1+percentage/100, that is, as a number greater than or equal to 1 and less than 2 (for the most part at least). As such, fractional representations of interest rates (usually given as a decimal number between 0 and 1 need to be converted before use. This function converts those decimal fractions into the right range.
percentage_to_rate(percent)
This library expects interest rates to be given to functions in the form of 1+percentage/100, that is, as a number greater than or equal to 1 and less than 2 (for the most part at least). As such, percentage representations of interest rates (usually given as a decimal number between 0 and 100 need to be converted before use. This function converts those percentages into the right range.
rate_to_fraction(rate)
Converts the interest rate format used by this library into a decimal fraction.
rate_to_percentage(rate)
Converts the interest rate format used by this library into a percentage.
compound_fc(ic,rate,nt)
Basic compound interest; Find final capital (fc) from initial (ic) at given rate and number of terms (nt). e.g. £50 at 3.4% for 20 years (interest added once per year) becomes compound_fc(50, 1.034, 20) which yields an answer of approximately £97.58
compound_ic_from_fc(fc,rate,nt)
Inverting basic compound interest; Find the initial capital (ic) given the final capital (fc), the interest rate and the number of terms. e.g. if after 20 years at 3.4%, we have £97.58, what was the initial investment. This becomes compound_ic_from_fc(97.58, 1.034, 20) which gives the answer of 49.9977, which is £50 when rounded up to the nearest penny.
compound_nt_from_fc(fc,ic,rate)
Inverting basic compound interest; Find the number of interest-adding terms if given the final capital (fc), initial capital (fc) and the interest rate.
compound_rate_from_fc(fc,ic,nt)
Inverting basic compound interest; Find the interest rate given the final capital (fc), the initial capital (ic) and the number of terms.
loan_paym(ic,rate,nt)
Loan amortisation; Determine the payment per term in order to pay off the initial capital (ic) of a loan given the interest rate and the number of terms (nt). Assumes that interest is always added before a payment is subtracted from the remaining capital, and that this is done only once per term. Given that terms are generally years, this is unusual, but not unheard of.
loan_apaym(*a__[],ic,rate,nt)
Loan amortisation; Uses the undocumented pass-by-reference feature of bc in order to create an array of term-by-term values showing remaining balance when the optimal payment is taken. To actually determine the optimal term payment, use the related loan_paym function which, other than the array reference, has the same parameter layout.
loan_paym_split(ic,rate,nt,spt)
Loan amortisation; Determine the payment per sub-term (sub-terms tend to be months within a year) in order to pay off the initial capital (ic) of a loan given the interest rate per term (terms tend to be years) and the number of terms (nt). Assumes that interest is always added before a payment is subtracted from the remaining capital. e.g. To calculate the monthly payment on a mortgage of £125,000 over 25 years at 5.66% APR with interest and payments being applied and taken monthly, use loan_paym_split(125000, 1.0566, 25, 12). This yields a suggested monthly payment of £768.98 when advantageously (for the borrower) rounding up to the next penny.
loan_apaym_split(*a__[],ic,rate,nt,spt)
Loan amortisation; Uses the undocumented pass-by-reference feature of bc in order to create an array of subterm-by-subterm values showing remaining balance when the optimal payment is taken. To actually determine the optimal subterm payment, use the related loan_paym_split function which, other than the array reference, has the same parameter layout.
loan_tpaym(ic,rate,nt)
Loan amortisation; Determine the total payment (tpaym) once the initial capital (ic) has been fully paid off when given the interest rate and the number of terms (nt). Assumes that interest is always added before a payment is subtracted from the remaining capital, and that this is done only once per term. Given that terms are generally years, this is unusual, but not unheard of.
loan_tpaym_split(ic,rate,nt,spt)
Loan amortisation; Determine the total payment (tpaym) once the initial capital (ic) has been fully paid off when given the interest rate per term, the number of terms (usually years) and number of subterms per term (usually 12 months). e.g. To calculate sum of all 300 monthly payments on a mortgage of £125,000 over 25 years at 5.66% APR with, use loan_tpaym_split(125000, 1.0566, 25, 12). This yields a value suggesting that over that time, the borrower will have to pay back a total of £230,692.14.
loan_ic_from_paym(paym,rate,nt)
Inverting loan amortisation; Given the term payment (paym), the interest rate and the number of terms, determine what the initial capital (ic) must have been.
loan_ic_from_paym_split(paym,rate,nt,spt)
Inverting loan amortisation; Given the sub-term payment (paym), the interest rate per term, and the number of terms and subterms per term (nt and spt), determine what the initial capital (ic) must have been.
loan_ic_from_tpaym(tpaym,rate,nt)
Inverting loan amortisation; Given the total payment over all terms (tpaym), the interest rate, and the number of terms (nt), determine what the initial capital (ic) must have been.
loan_ic_from_tpaym_split(tpaym,rate,nt,spt)
Inverting loan amortisation; Given the total payment over all sub-terms (tpaym), the interest rate, the number of terms and subterms per term (nt and spt), determine what the initial capital (ic) must have been.
loan_nt_from_paym(paym,ic,rate)
Inverting loan amortisation; Given the preferred payment per term (paym), the initial capital borrowed (ic) and the interest rate, determine how many terms are required to pay off the loan. Given that terms are usually years, this is somewhat unusual but terms can also be months if a monthly payment is given.
loan_nt_from_paym_split(paym,ic,rate,spt)
Inverting loan amortisation; Given the preferred payment per subterm (paym), the initial capital borrowed (ic), the interest rate and the number of subterms per term (spt), determine the number of full terms are required to pay off the loan.
loan_nt_from_tpaym(tpaym,ic,rate)
Inverting loan amortisation; Given the maximum preferred total payment over all terms (tpaym), the initial capital borrowed (ic) and the interest rate, determine how many terms are required to pay off the loan.
loan_nt_from_tpaym_split(tpaym,ic,rate,spt)
Inverting loan amortisation; Given the maximum preferred total payment over all subterms (tpaym), the initial capital borrowed (ic), the interest rate and the number of subterms per term (spt), determine how many full terms are required to pay off the loan.
loan_rate_from_paym(paym,ic,nt)
Inverting loan amortisation; Given the preferred payment per term (paym), the initial capital borrowed (ic) and the required number of terms (nt), determine the optimal interest rate which best fits these options.
loan_rate_from_paym_split(paym,ic,nt,spt)
Inverting loan amortisation; Given the preferred payment per subterm (paym), the initial capital borrowed (ic), the required number of terms (nt), and the number of subterms per term (spt) determine the optimal interest rate (for a term; usually a year) which best fits these options.
loan_rate_from_tpaym(tpaym,ic,nt)
Inverting loan amortisation; Given the preferred total payment over all terms (tpaym), the initial capital borrowed (ic) and the required number of terms (nt), determine the optimal interest rate which best fits these options.
loan_rate_from_tpaym_split(tpaym,ic,nt,spt)
Inverting loan amortisation; Given the preferred total payment over all subterms (paym), the initial capital borrowed (ic), the required number of terms (nt), and the number of subterms per term (spt) determine the optimal interest rate (for a term; usually a year) which best fits these options.
loan_spt_from_paym(paym,ic,rate,nt)
Inverting loan amortisation; Given the preferred payment per subterm (paym), the initial capital borrowed (ic), the interest rate and the number of terms (nt), determine the required number of subterms per term (spt) in order to correctly pay off the loan. This function is unlikely to see much use, except as a curiosity; Most loans are paid back monthly, meaning that given most actual loan data, this is likely to return an answer somewhere around 12.
loan_spt_from_tpaym(tpaym,ic,rate,nt)
Inverting loan amortisation; Given the preferred total payment over all subterms (tpaym), the initial capital borrowed (ic), the interest rate and the number of terms (nt), determine the required number of subterms per term (spt) in order to correctly pay off the loan. As above, this function is included as a curiosity, and will return an answer somewhere around 12 (months per year) when given actual loan data.
saving_fc(ic,paym,rate,nt)
Savings with compound interest; Given a starting amount, that is, some initial capital (ic), as well as a term-wise (terms here are usually years) single saving payment (paym) the financial institution's interest rate and a number of payment terms (nt), determine the final capital (fc) after those terms are over.
saving_afc(*a__[],ic,paym,rate,nt)
Savings with compound interest; Uses the undocumented pass-by-reference feature of bc in order to create an array of term-by-term values showing current balance after interest and payment have been added. Is otherwise identical to the saving_fc() function.
saving_fc_split(ic,paym,rate,nt,spt)
Savings with compound interest; Given a starting amount, that is, some initial capital (ic), as well as a subterm-wise (subterms here are usually months) single saving payment (paym) the financial institution's interest rate per term (~yearly) the number of payment terms (nt), and the number of subterms per term (usually 12 if monthly) determine the final capital (fc) at the end of that time. e.g. A prudent saver opens a 3.2% savings account with a lump sum of £5,000, and then pays in £260 per month over ten years. To find the final capital amount, we use saving_fc_split(5000, 260, 1.032, 10, 12) which reveals that the final amount should be somewhere around £43,476.14. With no interest, this would be merely £5,000 + 10*12*£260 or £36,200.
saving_afc_split(*a__[],ic,paym,rate,nt,spt)
Savings with compound interest; Uses the undocumented pass-by-reference feature of bc in order to create an array of subterm-by-subterm values showing current balance after interest and payment have been added. Is otherwise identical to the saving_fc_split() function.
saving_ic_from_fc(fc,paym,rate,nt)
Inverse calculations for savings with compound interest; Given the final capital (fc) required, the preferred term-wise payment (paym), the financial institution's interest rate and the preferred number of terms that the savings are to occur over (nt), determine what initial capital (ic) would be required in order to make it a reality.
saving_ic_from_fc_split(fc,paym,rate,nt,spt)
Inverse calculations for savings with compound interest; Given the final capital (fc) required, the preferred subterm-wise payment (paym), the financial institution's interest rate, the preferred number of terms that the savings are to occur over (nt), and the number of subterms per term (spt) determine what initial capital (ic) would be required in order to make it a reality. e.g. in a previous example, an investor started with initial capital of £5,000; If they wished to have exactly £50,000 at the end of the ten years, how much should they have started with? This is answered with saving_ic_from_fc_split(50000,260,1.032,10,12) yielding an answer of £9761.11 when rounded up to the next penny.
saving_nt_from_fc(fc,ic,paym,rate)
Inverse calculations for savings with compound interest; Given the required final capital (fc), and the initial capital (ic) as well as the preferred term-wise payment and the financial institution's interest rate, determine how many terms (nt), that is, how long to continue investing, until the requirement is met.
saving_nt_from_fc_split(fc,ic,paym,rate,spt)
Inverse calculations for savings with compound interest; Given the required final capital (fc), and the initial capital (ic) as well as the preferred term-wise payment, the financial institution's interest rate and the number of subterms per term (spt) determine how many terms (nt), that is, how long to continue investing, until the requirement is met.
saving_paym_from_fc(fc,ic,rate,nt)
Inverse calculations for savings with compound interest; Given the required final capital (fc), and the initial capital (ic) as well as the financial institution's interest rate and the number of terms of investment (nt), determine the optimal term-wise payment (paym) in order to meet the target. This is unusual as it assumes one lump sum per term, and terms here are usually considered to be years.
saving_paym_from_fc_split(fc,ic,rate,nt,spt)
Inverse calculations for savings with compound interest; Given the required final capital (fc), and the initial capital (ic) as well as the financial institution's interest rate, the number of terms of investment (nt), and the number of subterms per term (spt) determine the optimal subterm-wise (usually monthly) payment in order to meet the target.
saving_rate_from_fc(fc,ic,paym,nt)
Inverse calculations for savings with compound interest; Given the required final capital (fc), and the initial capital (ic) as well as the preferred term-wise payment (paym) and number of terms of investment (nt), determine what interest rate would be required from the financial institution in order to reach the goal. [A pipe dream for sure, as no-one gets to set their own interest rate!]
saving_rate_from_fc_split(fc,ic,paym,nt,spt)
Inverse calculations for savings with compound interest; Given the required final capital (fc), and the initial capital (ic) as well as the preferred term-wise payment (paym),number of terms of investment (nt), and number of subterms per term (usually 12 for months per year) determine what term-wise interest rate would be required from the financial institution in order to reach the goal.
saving_spt_from_fc(fc,ic,paym,rate,nt)
Inverse calculations for savings with compound interest; Another curiosity function included for the sake of completeness. Given preferred final capital (fc), initial capital (ic), preferred subterm-wise payment, an interest rate and a number of terms of investment (nt) determine the number of hypothetical payment and interest increment steps, that is, subterms per term (spt) required in order to reach the goal. Most ordinary financial data will cause this function to return 12 (months per year), but pathological cases may cause strange errors or impossible answers.
logic.bc
A large suite of functions to perform bitwise functions such as AND, OR, NOT and XOR. Uses twos complement for negative numbers, unlike previous versions of this file, which had no support at all.

Some of the functions here will use the global bitwidth variable, which itself is initialised as part of this file, to emulate byte/word sizes found in most computers. If this variable is set to zero, an infinite bitwidth is assumed.

Many functions will display a warning if there is suspicion that a secondary floating point representation of a number has been generated, e.g. 1.1111... is an SFPR of 10.0000...; These warnings can be disabled by setting the global variable sfpr_warn to 0 (default is 1).

bitwidth(x)
This function determines the minimal bitwidth needed to contain the value of x. Effectively an integer logarithm function.
bw_mult_(sc)
Internal function: Used along with internal global variables bw_mult_ml_ and bw_mult_sc_ to help manage the floating point bitwise functions.
checkbitwidth_()
Internal function: Used to check, before use, that the global bitwidth variable has not been set to an invalid value.
and(x,y)
Performs a bitwise logical AND of x and y. Uses base 4 internally for faster calculation.
andf(x,y)
As above but includes any floating point portion which may be present.
bitrev(x)
Reverse the bits in x. Uses bitwidth if it is nonzero.
graycode(x)
Convert x into its Gray code equivalent
graycodef(x)
As above but includes any floating point portion which may be present. NB: Since floating point allows carries of bits over to fractional bit positions, this function will not necessarily return the same answer as the above, being greater by 0.5 in those cases
hamming(x,y)
Find the binary Hamming distance between two numbers.
inverse_graycode(x)
Converts Gray encoded x back into its original bit pattern
inverse_graycodef(x)
Floating point inverse Gray code. All the caveats of the graycodef() function apply.
is_sfpr_(x)
Internal function to determine whether x is a secondary floating point representation (see above).
is_any_sfpr3_(x,y,z)
Internal function to determine whether one or more of x, y or z is a secondary floating point representation (as above).
not(x)
Perform a bitwise logical NOT of x. Since these functions use twos complement, this function returns -1-x which has an exactly flipped bit representation in 2C.
or(x,y)
Perform a bitwise logical OR of x and y. Uses base 4 internally for faster calculation.
orf(x,y)
As above but includes any floating point portion which may be present.
orm(x,y)
'Multiplies' x and y in a no-carry, bitwise manner using logical OR in place of addition.
ormf(x,y)
As above but includes any floating point portion which may be present.
resign(x)
Despite an apparently pessimistic name, this actually RE-applies a SIGN to x, with the assumption that the current bitwidth is valid. e.g. if bitwidth is 8, resign(254) is -2. C programmers will recognise this as effectively casting an unsigned value into a signed variable of the same size.
rol(x,n)
Roll Left: Familiar to assembly programmers, this shifts x left by n places within the current bitwidth and adds the carried left hand bits back on the right. e.g. 10010011 rolled left by 3 is 10011100 assuming a bitwidth of 8.
ror(x,n)
Roll Right: As above but shifts to the right, placing lost right hand bits back on the left. May well complain if bitwidth is 0 (i.e. implied infinite), as right hand bits would have to be placed in infinite positions.
sfpr_warn_msg_()
Internal function to display the aforementioned warning about SFPRs.
shl(x,n)
Shift Left: Shifts the bits in x left by n places. Bits carried from the left hand side are lost if x cannot be kept within the current bitwidth.
shr(x,n)
Shift Right: Shifts the bits in x right by n places. Bits from the right hand side are lost
ungraylike1(x)
Self-inverse binary permutation of x. At first glance resembles a Graycode style transformation, but this is not the case, hence "ungraylike". Reverses the order of all bits after the first.
ungraylike2(x)
Another self-inverse binary permutation of x. Also resembles a Graycode style transformation, but is also not the case. Reverses the order of and flips all bits after the first.
unsign(x)
Interpret a negative number as a positive number within the current bitwidth. Again, C programmers will recognise this as a cast from signed to unsigned.
xor(x,y)
Perform a bitwise logical XOR (EXCLUSIVE OR) of x and y. Uses base 4 internally for faster calculation.
xorf(x,y)
As above but includes any floating point portion which may be present.
xorm(x,y)
'Multiplies' x and y in a no-carry, bitwise manner using logical XOR in place of addition.
xormf(x,y)
As above but includes any floating point portion which may be present.
logic_andm.bc
Various attempts to create bitwise AND 'multiplication' functions that do not result in zero
andm(x,y)
'Multiplies' x and y in a no-carry, bitwise manner using logical AND. This function would return zero all the time but has been tweaked to return values where possible by starting out with a substrate of 1s.
andmf(x,y)
As above but includes any floating point portion which may be present.
x1andm(x,y)
One method of combining the equivalences between AND, OR & XOR alongside their multiplicative equivalents in order to emulate AND-multiplication
x1andmf(x,y)
As above but includes any floating point portion which may be present.
x2andm(x,y)
Second method of combining the equivalences between AND, OR & XOR alongside their multiplicative equivalents in order to emulate AND-multiplication
x2andmf(x,y)
As above but includes any floating point portion which may be present.
logic_inverse.bc
New for February 2013
Calculation of possible inputs to bitwise AND and OR functions in order to achieve a specific output, i.e. find inverses, where possible, to bitwise functions.
dand_sor_(which, z,y ,n)
Internal function; The engine which finds all inverse solutions in this library.
dand(z,y ,n)
"Division AND" / "De-AND"; Attempts to find solutions to bitwise z = and(x,y) for x. Often, there are multiple solutions or even none at all and so the n parameter determines the kind of return value to be given by the function. This can either be a coded representation of all possible solutions; one of the solutions, should one exist; or even an error code of -1 if there are no solutions. A warning will also be issued if there are no solutions when a single solution is requested and the dand_sor_warn_ variable is non-zero. This is enabled by default. Values for n are:
n = -4
Requests the return of a base 4 (quartal) codification of all solution bits. For the quartal digits in the return value:

0 represents the fact that no solution bit is possible for the respective binary position in a potential solution. Any zero quartal digits in this codification therefore mean that no solution is possible for the other parameters.

1 represents that a zero bit and only a zero bit is possible at the respective binary position in a potential solution.

2 represents that a one bit and only a one bit is possible at the respective binary position in a potential solution. Note that in this representation, zero and one are not represented by their own values, despite there otherwise being some logic to the assignation of the quartal digits.

3 represents a don't-care state in that either a zero or a one bit is suitable for the respective binary position in a potential solution. For every 3 in this representation, the potential number of solutions doubles.

n = -3
Requests the return of a base 4 (quartal) codification of all solution bits. For the quartal digits in the return value:

0 represents that a zero bit and only a zero bit is possible at the respective binary position in a potential solution.

1 represents that a one bit and only a one bit is possible at the respective binary position in a potential solution. Note that in this representation, zero and one are represented by their own values, making any direct quartal output easier to read.

2 represents a don't-care state in that either a zero or a one bit is suitable for the respective binary position in a potential solution. For every 2 in this representation, the potential number of solutions doubles.

3 represents the fact that no solution bit is possible for the respective binary position in a potential solution. Any '3' quartal digits in this codification therefore mean that no solution is possible for the other parameters.

This quartal representation choice is the default in this library, not least because the mnemonic is simple: 0 is 0, 1 is 1, and 2 doubles the solution count.

n = -2
Requests the return of the number of possible solutions. Doing this will produce no warnings as zero solutions is a valid and expected answer in some cases.
n = -1
Requests the return of the number of don't-care states, which is log2 of the number of solutions. Also does not warn if there are no solutions, but will cause the return of the more sensible -1, rather than negative infinity, when no solutions exist.
n ≥ 0
Requests the return of the nth solution (Counting from 0 and modulo the number of solutions. This will be a power of two if any exist.) for the other parameters. If there are no solutions, a warning will result (assuming the dand_sor_warn_ variable is non-zero) and -1 will be returned instead of a solution.

Note that setting n to 0 guarantees that a solution is returned should one exist.

sor(z,y ,n)
"SubtractOR" / "SubtORct"; Attempts to find solutions to bitwise z = or(x,y) for x. Often, there are multiple solutions or even none at all and so the n parameter determines the kind of return value to be given by the function. See the definition of dand() for an explanation of the values n can take.
striped_dand(z,y ,n)
Striped "Division AND" / "De-AND"; Attempts to find solutions to bitwise z = striped_and(x,y) for x. Often, there are multiple solutions or even none at all and so the n parameter determines the kind of return value to be given by the function. See the definition of dand() for an explanation of the values n can take.

For an explanation and definition of "striped AND", see this site's logic_striping.bc section.

striped_sor(z,y ,n)
Striped "SubtractOR" / "SubtORct"; Attempts to find solutions to bitwise z = striped_or(x,y) for x. Often, there are multiple solutions or even none at all and so the n parameter determines the kind of return value to be given by the function. See the definition of dand() for an explanation of the values n can take.

For an explanation and definition of "striped OR", see this site's logic_striping.bc section.

print_01dx_(x)
Internal output driver for interpreting and then printing the default n = -3 quartal solutions pattern.
dand_print(z,y)
Prints a specially formatted quartal string to aid in identifying, by sight, any and all possible solutions to bitwise z = and(x,y) for x. Output contains "0" for "position must be 0", "1" for "position must be 1", "d" for "don't care / 0 or 1 is fine" and "X" for "impossible / fail bit". Any "d"s in an output indicate multiple solutions. Any "X"es indicate that no solution is possible.
sor_print(z,y)
Prints a specially formatted quartal string to aid in identifying, by sight, any and all possible solutions to bitwise z = or(x,y) for x.
striped_dand_print(z,y)
Prints a specially formatted quartal string to aid in identifying, by sight, any and all possible solutions to bitwise z = striped_and(x,y) for x.
striped_sor_print(z,y)
Prints a specially formatted quartal string to aid in identifying, by sight, any and all possible solutions to bitwise z = striped_or(x,y) for x.
logic_otherbase.bc
An abortive attempt to scale up bitwise functions to bases other than binary; The well defined and more complete set of binary bitwise functions can be found in logic.bc, although all functions herein will return correct results when base 2 is specified.
asym_parity(base,x,y)
A possible XOR extension; Inverts (i.e. subtracts from base-1) those digits in x which are flagged as to be flipped by the digits in y. In even bases, this amounts to "invert digit in x if digit in y is odd." Is asymmetric in that asym_parity(base, x,y) almost never equals asym_parity(base,y,x).
asym_mixor(base, x,y)
A second, similarly asymmetric possible extension of bitwise XOR into other number bases; A blend of no_borrow_diff() and no_carry_add().
base_graycode(base,x), inverse_base_graycode(base,x)
An extension to Graycode into other number bases, preserving the requirement that in the sequence of Gray codes, no two adjacent codes differ by more than one digit.
base_hamming(base,x,y)
As with binary, this is a difference-counting function, returning the total number of differences between two numbers in a given base. e.g. the number of differences between 1234 and 9284 in base ten is 2, since they differ in the first and third positions (when read left-to-right). For an alternative to this function which sums the actual difference, which would be thirteen in this example, see the digit_distance() function in digits.bc.
digitwise_diff(base, x,y)
Another possible extension of XOR into other number bases. Akin to, but not to be confused with the above, as this returns a number rather than a sum of differences.
digitwise_sdiff(base, x,y)
A variant of the above, so XOR again. In modulo arithmetic, there are two solutions to finding the difference between two numbers. e.g. |1-7| modulo ten could be 6 or 4. This function chooses the lower of the two.
digitwise_max(base, x,y)
A logical extension of bitwise OR.
digitwise_min(base, x,y)
A logical extension of bitwise AND.
digitwise_modmult(base, x,y)
Another logical extension of bitwise AND.
digitwise_tlumdom(base, x,y)
Another logical extension of bitwise OR. Note that tlumdom is modmult backwards; the underlying algorithms are related.
digitwise_xor_(which, base, x,y)
Internal function; Contains the main engine for all of the XOR related functions in this file.
no_borrow_diff(base, x,y)
Yet another XOR-like extension. Takes the difference of digits of x and y, modulo the base, without borrowing from the left.
no_carry_add(base, x,y)
A sixth XOR-like extension. Adds digits of x and y, modulo the base, without carrying any overflow leftwards.
logic_striping.bc
A family of functions related to the bitwise functions which may be useful for encryption and hashing. Then again they might not. See the text documentation for more technical information.

These were separated from the main logic.bc due to being of questionable worth, but were given their own file as they are still interesting functions. The latter file is required for these functions to work correctly.

stripe_(b,x,y)
Internal function: Engine for striped_and and striped_or
striped_and(x,y)
Performs a bitwise logical 'STRIPED AND' of x and y
striped_andf(x,y)
As above but includes any floating point portion which may be present.
striped_andm(x,y)
'Multiplies' x and y in a no-carry, bitwise manner using logical 'STRIPED AND' in place of addition.
striped_andmf(x,y)
As above but includes any floating point portion which may be present.
striped_or(x,y)
Performs a bitwise logical 'STRIPED OR' of x and y
striped_orf(x,y)
As above but includes any floating point portion which may be present.
striped_orm(x,y)
'Multiplies' x and y in a no-carry, bitwise manner using logical 'STRIPED OR' in place of addition.
striped_ormf(x,y)
As above but includes any floating point portion which may be present.
genstripe(override,repeat,x,y)
A generalisation of the concept behind 'STRIPED AND' and 'STRIPED OR' creating an infinite class of related bitwise functions.
genstripef(o,r,x,y)
As above but includes any floating point portion which may be present.
genstripem(override,repeat,x,y)
'Multiplication' as seen many times before here.
genstripemf(o,r,x,y)
As above but includes any floating point déjà vu which may be present.
logic_striping_meta.bc
New for February 2013
A library for exploring the mathematics of the striping pattern numbers used in logic_striping.bc which allows integers to be interpreted in a completely different way to either integer or modular arithmetic. See the end of the text documentation for more technical information.
simplify_stripe_pattern(x)
Converts x into the minimal stripe pattern value. e.g. if x, interpreted in binary, is of form 1[pattern][pattern]... where the same pattern (bit length unimportant) repeats to the end of the number, then this function will return the number whose binary representation is 1[pattern]. e.g. when given the number 26dec = 11010bin, this function detects and removes the repeat of '10' at the end and returns 110bin = 6dec, likewise, 731dec = 1011011011bin → 1011bin = 11dec. This is somewhat similar to finding the smallest congruence in modular arithmetic, or a prime given a power of that prime.
next_match_stripe_pattern(x)
Generates the next stripe pattern in a sequence by appending one iteration of any repeating bit pattern found after the leading 1. If the bit pattern does not repeat, the number returned will have two copies of the whole bit string part of its pattern interpretation. e.g. 1011011 → 1011011011 and 1010011010 → 1010011010010011010. This is somewhat similar to finding the next largest congruence in modular arithmetic, or the next largest power of a prime given particular power of that prime.
rep_stripe_pattern(x,p)
Takes the entire bit pattern (repeating or otherwise) after the leading 1 in the binary representation of x and repeats it p times. e.g. 1010010 → 1010010010010. Note that no simplification of the bit pattern is performed. This has similarities to raising a number to a particular power or multiplying unary numbers.

Negative values of p will invert the bit pattern after the leading 1 before repeating the pattern the specified number of times. Indeed all functions here will accept a negative value for their pattern parameter, interpreting the pattern as its bit-flipped equivalent. e.g. -1010010 ↔ 1101101. This means that there are technically no 'negative' patterns as they all have a positive interpretation, similar to how negative numbers are merely an interpretation of a residue in modular arithmetic.

repsof_stripe_pattern(x)
Determines the maximal number of repeats of the bit pattern after the leading 1 in the binary representation of x. e.g. 1010010010010 has four copies of the substring '010', and so this function would return 4 given this number. Note that this could also have been interpreted as two repeats of '010010', but this is not a maximal solution. This is similar to determining the power of a prime when given only the value of that prime power.
stripe_pattern_to_1c(x)
Converts a stripe pattern into a 1s complement integer such that the integer is simply the stripe pattern without the leading 1. Since this would cause ambiguity for those patterns which begin with a zero (01 = 001 = 0001 etc.), all patterns with a leading zero are bit flipped and negated. e.g. pattern value 11101 becomes 1101, but pattern value 10101 becomes -1010. Note that this is very similar to the procedure for negative values as described within the definition for rep_stripe_pattern().
stripe_pattern_to_2c(x)
As above, but negative values are calculated as if they are 2s complement, meaning that when read as a binary value they are one away from the value of the pattern they represent. e.g. 1011011011 → -100100101
stripe_pattern_from_1c(x)
Inverse of stripe_pattern_to_1c(); Creates a standard stripe pattern integer from a 1s complement representation.
stripe_pattern_from_2c(x)
Inverse of stripe_pattern_to_2c(); Creates a standard stripe pattern integer from a 2s complement representation.
mul_stripe_patterns(x,y)
A method for combining / multiplying two patterns. The result is equal to the first pattern repeated either as-is, or bit-flipped, designated by the bits of the second. This multiplication is not symmetric except if the parameter patterns are powers of two or one less than a power of two. Indeed this multiplication implies a direct connection between these numbers and the positive and negative integers. See comments in code. The number of bits in the result is equal to the product of the number of bits in the parameters.
div1_stripe_patterns(z,y)
Since the above multiplication is not symmetric, there are two division functions. Given the result of such a multiplication and the right hand side value of the multiplication, this reconstructs the left hand (or 1st) value.
div2_stripe_patterns(z,x)
Since the above multiplication is not symmetric, there are two division functions. Given the result of such a multiplication and the left hand side value of the multiplication, this reconstructs the right hand (or 2nd) value.
sqrt_stripe_pattern(z)
Given the result of one of the above multiplications, try to find whether it could have been the same pattern multiplied by itself.
mix_stripe_patterns(x,y)
An alternative and symmetric multiplication of patterns using NXOR. The length of the resulting pattern is the lowest common multiple (not the product) of the lengths of the two parameters. Multiplicative-like structure is conserved as is the sign of the result when considering the previously mentioned mapping of powers of two and one less than powers of two to negative and positive integers.
unmix_stripe_patterns(z,x)
The inverse of the above up to pattern uniqueness in the sense of representing the same repeating bit string. Given a result of the mix-multiplication and one of the parameters, reconstruct the other.
modsum_stripe_patterns(x,y)
Another method for combining patterns using modular addition, modulo 2lcm of pattern lengths.
unmodsum_stripe_patterns(z,x)
The inverse of the above up to pattern uniqueness in the sense of representing the same repeating bit string. Given a result of the modular sum and one of the parameters, reconstruct the other.
cat_stripe_patterns(x,y)
The most basic way of combining patterns, akin to addition; Create a new pattern by appending one to another. (Cat = catenation). Has the unfortunate side effect of not properly preserving addition for the aforementioned powers of two and one less than powers of two. e.g. Adding the pattern akin to -3 to the pattern akin to 2 will not result in the pattern akin to -1.
decat_stripe_patterns(z,y)
An extended inverse of the above. Given a pattern and a supposed right hand parameter from a prior catenation, chop off the matching segment of the catenation, or else catenate the bit-flipped version of the right hand parameter. See notes in code.
decat2_stripe_patterns(z,y)
A more intelligent version of the above. Removes the largest matching part of the right hand side of the left parameter relative to the left hand side of the right parameter's pattern. Any non-matching remainder of the right hand parameter is then bit flipped and appended to the result as in the previous decatenation. See notes in code.
undecat2_stripe_patterns(z,y)
The pièce-de-résistance; This function preserves the addition of the power-of-two / one-less-than-power-of-two patterns when interpreted as negative and positive integers. Removes the largest matching-when-bit-flipped portion of the right hand side of the left parameter relative to the left hand side of the pattern of the right parameter, before appending the non-bit-flipped unmatching portion of the right parameter. See notes in code.
melancholy.bc and melancholy_b.bc
Twinned for February 2013
Twin suites of functions for investigating the two kinds of Melancholy Numbers and the iterations which lead to them. These are a discovery undoubtedly made by many people, myself included.

Your humble author is guilty of coining the term "melancholy" with the intention of comparing these with the unhappy numbers. See here for a Usenet discussion on the first and original kind, and also many notes stored within the bc code in both files.

is_melancholy(x)
is_melancholyb(x)
Returns 1 (true) if x is melancholy, i.e. the iteration does not reach 0, and 0 otherwise.
melancholy_chainlength(x)
melancholyb_chainlength(x)
Returns the number of iterations required to reach 0 for the given x. Will return a negative number of iterations if x is melancholy, which is the number of iterations required before the algorithm "realises" that it will not reach 0.
melancholy_lastsqrt(x)
melancholyb_lastsqrt(x)
The nature of the melancholy algorithm is such that a square number immediately leads to 0, meaning that x is not melancholy. This function gives the square root of that final perfect square.
melancholy_loopsize(x)
melancholyb_loopsize(x)
Returns the number of iterations in the loop encountered by numbers that are melancholy.
melancholy_max(x)
melancholyb_max(x)
During determining whether the iteration of x will eventually reach 0, the iterates rise and fall much like those of the Collatz conjecture. This returns the largest number encountered before reaching zero, or else prior to detecting a loop.
melancholy_print(x)
melancholyb_print(x)
Shows all the iterations down to 0 if the number is not melancholy, or else stops once a loop is detected.
melancholy_root(x)
melancholyb_root(x)
Shows 0 if the number is not melancholy, or else the smallest number in the loop that a melancholy number becomes trapped within.
is_melancholy_sg(x)
is_melancholyb_sg(x)
sg = "set globals": This function is all of the above rolled into one, and will set global variables by the names of the above functions (e.g. melancholyb_max) for the parameters given. f the global variable melancholy_print is set to 1, then this function will also behave as the melancholy_print() function and display the value of the iterations. Set to 0 to turn the feature off again. Like is_melancholy(), returns 1 if x is melancholy and 0 otherwise.
misc_235.bc
Find the sum of powers of 2, 3 and 5 that are closest to a number
print235(x)
The only function in this odd little file. Does exactly as described, although the implementation isn't perfect and it sometimes misses a nearer answer to the one given if there isn't an exact solution.
misc_ack.bc
Calculate the hyper-exponential Ackermann function; All but useless given that bc can't cope with such huge numbers!
ack(x,y)
Tries to calculate the Ackermann function
ackz(x,y)
As above, but works better (when it works at all) for floating point values
misc_anglepow.bc
A semi linear version of exponentiation and logarithm related to the number sequence: 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 2000, etc. [See the OEIS's A037124]
anglepow10(x)
Calculates the x'th entry in the above sequence
anglelog10(x)
Calculates the inverse of the above
anglepow(b,x)
An extension of anglepow10 to all number bases b
anglelog(b,x)
Calculates the inverse of the above
misc_perfectpow.bc
Find the nearest perfect power to a number
nearest_perfect_power(x)
Finds the nearest perfect power to x, e.g. the nearest perfect power to 31 is 32 which is 25. In the case of a tie, chooses the lower option, e.g. 26 is midway between 52 = 25 and 33 = 27 and so 25 is returned.
nearest_perfect_power_a(*a__[],x)
As above, but uses the undocumented pass-by-reference for arrays to store extra information into the given array. Array element 0 is set to the perfect power itself, elements 1 and 2 are the root of the power and it's radix, and the third element is the sign of x. e.g. x = -124 would cause the array to contain {125,5,3,-1} because 125 = 53 and the input was negative.
misc_srr.bc
The Sum of Repeated Roots function
The mathematical formula for this interesting function is:
SRR(x) =
∞
Σ
n = 1
 x2-n - 1

Or as simply as possible: srr(x) = sum[n=1..oo] x^(2^(-n))-1
srr(x)
The function as described above.
srr_n(n,x)
A generalisation of srr, which uses the nth root rather than the square root. srr_n(2,x) is equivalent to srr(x)
orialc.bc
Extensions and unusual variants of factorial-related functions. Requires all three of funcs.bc, factorial.bc and primes.bc in order for all functions to work.
lcmultorialc(x)
Simple and logical extension to the lcmultorial() as found in factorial.bc. Simply multiplies by a fractional power of the next LCMultorial.
primorialc(x)
Dummy function which suggests the use of any of the other primorialc alternatives.
primorialc_fact(x)
An attempt to extend the primorial function over the reals by using the factorial() function as a substrate.
primorialc_nextp(x)
A very simple extension to the primorial which multiplies by a fractional power of the next prime when given a non-prime argument.
primorialc_self(x)
Another simple extension to the primorial; Multiplies by a fractional power of the argument when it is not prime.
primorialc_backstep(x)
A more complicated combination of the above two functions. Calculates the next highest primorial and then uses the parameter to determine which fractional power of parameter to divide by to find a sensible intermediate value.
primorialc_accident(x)
An incorrect version of the above which has a serendipidous bug leading not only to sensible values for the extended primorial, but also has some instances for non-prime integer argument where the result is actually rational. Examples may be found in the source code.
submodorial(x)
Submodulus superprimorial / Product of remainders / Product of x mod k for 2 ≤ k < x See A180491 in the OEIS.
submodorialg(n,x)
Generalised version of the above; Product of n+(x mod k) for 2 ≤ k < x
lcmsubmodorial(x)
LCM of remainders / LCM of x mod k for 2 ≤ k < x. LCM equivalent of submodorial.
lcmsubmodorialg(n,x)
Generalised version of the above; LCM of n+(x mod k) for 2 ≤ k < x. LCM equivalent of submodorialg.
output_formatting.bc
Updated for June 2012
Powerful formatting tools for GNU bc. Most functions here should have their return value assigned to a variable, and unless otherwise specified here, will return the value of the number they were asked to print.
comma_(x,gp)
Internal function: Engine for commaprint.
commaprint(x,g)
Print x with groups of digits of size g, separated by commas. e.g. commaprint(1020304,3) prints 1,020,304
intprint(w, n)
Print integer n within a field width of w characters.
letterl__(a), letteru__(a), letter_(mode, a)
Internal functions which print specific letters or digits given particular values of a. Lower case; upper case; mode sets alphanumeric or alpha-only.
newline()
Prints a newline
printbase(x)
By default bc will use decimal groups for 'digits' when outputting numbers with obase set above 16. When obase ≤ 16, letters are used as digits. This function outputs using letters right up to base 36 which uses Z as the one-less-than-base digit. Uses obase, so no base need be specified as a parameter.
Lowercase letters can be specified by setting the global variable output_lcase_ to 1. The default is 0, since bc usually prints its letter digits in uppercase.
printbase_letters(x)
Similar to the above, but uses the underscore symbol '_' for zero and the letters A to Z for the digits 1 onwards. This means that it will only work correctly for bases up to 27, at which point it defers to the above. Again, setting output_lcase_ to 1 will enforce lowercase.
printbijective(bbase,x)
Display x in the bijective base specified by bbase, which may be positive or negative. A special format has been devised to show - albeit invalidly - the fractional part of x. The value of zero is represented by a single dot, which doubles as the invalidating, bijection-breaking, basimal point. Note that printbijective(), unlike printbase(), requires the specification of the output base. bc's own obase, used by printbase(), cannot support base 1 (unary) whereas this function is perfectly capable, therefore the base must be specified.
printbijective_letters(bbase,x)
Similar to the above, the letters A to Z for the digits 1 onwards. This means that it will only work correctly for bases up to 26.
printsbase(base,x)
An all-in-one function that chooses which of the above four functions to call based on various global variables and settings. The first of these is the global bijective variable which is also found in digits.bc; When set to 1, it indicates bijective numeration output. The second is the global variable printsbase_letters_ which if set to 1, indicates that the output should be in letter form only.
printsobase(x)
As above, with the exception that the function will use bc's own obase global variable to specify the output base.
printfactorialbase(x)
Show the representation of x in the factoradic or factorial number base. Uses global variable pfactb_zero_ to specify whether the zeros for the 1! and 1/1! digit places are displayed. This is set to 0 by default and so these zeros are not displayed. Spaces are placed between digit places and individual digits are displayed in the current obase.
printnegabase_(base,x)
Internal function; Used by the above functions with a base or bbase parameter to allow them to support printing of numbers in negative bases.
printdms(x)
Treat x as a number of hours and print it in hh:mm:ss format.
printff(width, precision, n)
A function styled after the C syntax: printf("%*.*f",width,precision,n); There are some minor differences however:
  • If width is negative, n is aligned to the left hand side of the field
  • If precision is a non-integer, and alignment is set to the right hand side (i.e. width is positive) leading zeroes will fill the field prior to n.
  • If precision is set to zero, no decimal point nor fractional part of n will be displayed.
  • If precision is 0.0 (bc can tell the difference between 0 and 0.0), the above two features are combined.
printfe(width, precision, n)
A function styled after the C syntax: printf("%*.*e",width,precision,n); In addition to having the same features as for printff, above:
  • If precision is negative, the exponent will be set in engineering notation, i.e. will always be a multiple of 3 (assuming decimal output) It will, for example, choose to print 123.4e+00 rather than 1.234e+02, or 56.124e+06 rather than 5.6124e+07. The exponent multiple is calculated from the current obase, and so is not always 3.
printfrac(improper, maxdenom, x)
Prints x as the most accurate fraction possible under the restraint of a maximum denominator. Can choose improper fraction style if required. Will always choose a/b (proper fraction )style for fractions less than one. Output can be copy/pasted back into bc as valid syntax. Returns the value of the fraction printed and not the original value of x.
printsft(a,b)
Prints a and b as an improper fraction in Smallest Fractional Terms. Requires the GCD function from funcs.bc
printspc(n)
Prints the specified number of spaces.
printtabs(n)
Prints the specified number of tabs.
trunc(x)
Returns x with all trailing zeroes truncated, working around bc's habit of keeping them stored in the variable. Will also round up trailing nines (or "base-minus-one"s in other bases). Of course this latter case means that the return value is not always guaranteed to be equal to x, but may well have rounded x to the "correct" value.
output_graph.bc
A rudimentary console-based graphics package. Uses a global array called screen[] to store a very simple 'bitmap' made up of characters. The x and y dimensions are set with global variables screen_x and screen_y.
or_(x,y)
Internal function: Bitwise OR. See logic.bc for the fully fledged version of this function.
screen_clear()
Blanks out the global screen[] array.
screen_axes(xx,yy)
Draws axes with an origin at coordinates (xx,yy) into screen[] Returns 0 if (xx,yy) is out of bounds and 1 otherwise
screen_plot(x,y,c)
Put the character specified by c into screen[] the coordinates (x,y). If c is negative, an attempt is made to combine any character already at (x,y) with the character specified by -c. e.g. screen_axes(xx,yy) uses this feature to combine the y-axis "|" character with the x-axis "-" character to form a "+" sign at the origin. Character values have been chosen for c so that reasonable combinations will form using negative c. Returns 0 if (x,y) is out of bounds and 1 otherwise.
screen_print()
Actually display the intended interpretation of the contents of screen[] onto the console.
screen_printchar_(c)
Internal function: Prints a character specified by c.
output_roman.bc
A one-function library for displaying numbers in Roman style.
printroman(n)
Outputs n to the console in Roman numerals. Uses N (nulla) for zero, bracketed notation for thousands, millions, etc, and fractions are given in duodecimal unciae notation ("S::.")
Can use the output_lcase_ global variable, also found in the main output_formatting.bc, to specify lowercase numerals. Set it to 1 to enable this mode. Default is 0.
primes.bc
Updated for February 2013
A mostly naive implementation of prime number handling, with one exception: Contains a relatively powerful primality checker.

Many functions can be given a speed boost by using one of the prime databases available on this site and running the fillprimearray() function in primes_db_code.bc. This will fill the global prime[] array which the other functions can use as a reference but all functions can operate without that add-on. A word of caution however, regardless of whether a prime database is loaded: The main factorisation engine is a glorified Eratosthenes' sieve and so any functions here which rely on that engine may take a while to run.

To offset this somewhat, the last factorisation is always stored in global array factorpow[] so follow-up calculations on the same number do not need to perform the same factorisation again.

arithmetic_derivative(x)
Calculate the so-called Arithmetic derivative of x, using new internal factorisation storage routines.
count_divisors(x)
Also known as the sigma or sigma-zero function (see below), returns the total number of divisors of x, including 1 and x itself.
count_factors(x)
Returns the number of prime factors of x, counting primes more than once if necessary.
count_unique_factors(x)
Returns the number of unique prime factors of x, i.e. ignoring powers.
core(x)
Find the squarefree core of x, i.e. the product of all primes which have an odd power in the factorisation of x.
divisors_print(x)
Prints a partially ordered list of the divisors of x
divisors_sp_(*divs[],x,print_)
Internal function used by the immediately above and below.
divisors_store(*d[],x)
Stores a partially ordered list of the divisors of x into an array given as the first parameter. Preference is given to divisors ≤√x if the array is too small to hold all divisors; Larger divisors can be reconstructed from these. This uses the undocumented pass-by-reference feature of bc.
fac_print(x)
Formerly known as fac(x); Prints the prime factorisation of x
fac_sp_innerloop_()
fac_sp_(*fp[],x,print_)
fac_store_(*fp[],m,p,c,print_)
factorpow_set_(fp[])
Internal functions; These three are called by each other as well as the immediately above and below (print and store) functions
fac_store(*fp[],x)
Stores the prime factorisation of x into an array given as the first parameter. As before, this uses the undocumented pass-by-reference feature of bc. The array format is {prime_factor,power,prime_factor,power,...}. This function is also used internally by many other functions for calculations on factorisations.
has_freedom(x,y)
Determine whether x is y'th power-free. Returns 0 if the number is not power-free, but returns non-zero otherwise (not necessarily 1). When y = 1, defaults to returning whether x is prime. When y = 2, returns the Mobius function. For y > 2, returns a logical but non-standard extension to the Mobius function.
int_gcd(x,y)
Integer-only greatest common divisor function. This can also be found in funcs.bc.
int_modpow(x,y,m)
Quickly determine the value of xy mod m for x, y and m all positive integers.
int_modpow_(x,y,m)
Internal function, used by the above.
is_perfect_totient(x)
Returns whether a number is a so-called perfect totient number.
is_perrin_pseudoprime(p)
Determine whether p is a Perrin pseudoprime. Also performs a very basic Mersenne primality check.
is_practical(x)
Returns whether a number is a so-called practical number.
is_prime(x)
A powerful primality checker which combines the pseudoprime tests herein to determine whether x is extremely probably prime. If zero is returned, x is most definitely not prime. May take a while to return 1 in certain cases, but is almost certain to be correct if it does. Estimates suggest less than 1 in 101000 candidates would be misidentified.

Be aware that this function takes advantage of all pseudoprime tests in this file and so accuracy is susceptible to how the rabin_miller_maxtests_ global variable is set.

is_prime_td(x)
Determine whether p is prime solely using trial division. Uses the Perrin test as a shortcut, but uses trial division to check. May need years to run, but is guaranteed to return a right answer.
is_rabin_miller_pseudoprime(p)
Determine whether p is a Rabin-Miller pseudoprime . Performs multiple RM tests, to achieve high certainty of primality. The global variable rabin_miller_maxtests_, when set to a positive value, determines the maximum number of RM tests this function performs. This reduces time to execute at the expense of accuracy. By default this is set to 0, meaning run until almost certain the number is prime.
is_small_division_pseudoprime(x)
Confirms that x is indivisible by a significant number of small primes
largest_prime_factor(x)
Returns the largest prime factor of x; A full factorisation has to be performed, so this has more in common with other functions here than with its counterpart smallest_prime_factor().
largest_prime_power(x)
Returns whichever prime factor of x, when raised to its power from the full factorisation of x, is the largest.
primorial(n)
Returns the product of all primes ≤ n. See orialc.bc for a few possible extensions of the primorial to all numbers.
mobius(x)
The Mobius function: Returns 0 if x is not prime, or else -1 or 1 depending on the powers of the primes in x. Is an alias for has_freedom(x,2).
nextprime(n)
Returns the nearest prime > n; Relies on the is_prime() function when searching for candidates, so may take a while to find an answer.
nextprime_ifnotprime(n)
Returns the nearest prime ≥ n
nearestprime(n)
Returns the nearest prime to n, or n if n is prime; Relies on the is_prime(), prevprime(), and nextprime() functions when searching for candidates, so may take a while to find an answer.
prevprime(n)
Returns the nearest prime < n; Relies on the is_prime() function when searching for candidates, so may take a while to find an answer.
prevprime_ifnotprime(n)
Returns the nearest prime ≤ n
printfactorpow(fp[])
Prints the contents of the array interpreted as a prime factorisation. Expects the array to be of the form {prime_factor,power,prime_factor,power,...} as created by the fac_store() function.
rad(x)
Integer radical function; Returns the largest power-free number which is a divisor of x.
ramanujan_c(q,n)
Srinivasa Ramanujan's sum
sigma(n,x)
Sum of divisors function, with the extension of raising the divisors to the power n before summation. When n is 0, this is equivalent to the count_divisors() function. When n is 1, it is equivalent to the sum_of_divisors() function below.
smallest_prime_factor(x)
Returns the smallest prime factor of x; If this value is all that is required, this function is often much faster than performing a full factorisation as it stops once the factor is found.
squarepart(x)
Determine the square part of x, i.e. find the largest square number which is a divisor of x.
sum_of_divisors(x)
Simple sum of divisors function. Is an alias for sigma(1,x).
totient(x)
Euler totient function; Returns the number of integers less than or equal to x which are coprime to x.
totient_itercount(x)
Returns how many times the Euler totient function must be applied to x in order to reach 1.
totient_itersum(x)
Returns the sum of intermediate terms as the Euler totient function is repeatedly applied to x in order to reach 1.
primes_db.bc (Deprecated)
A database of primes in an array called prime[]. Contains all primes from the first to the 65,535th. This idea was independently discovered, but is identical to a technique found in the X-bc project.

To distinguish their project from this one, this file was generated by native, independent bc code and is "compressed" through use of hexadecimal.

NB: Even though this file is accurate, it is now deprecated and one of the other prime databases should be used instead.

* This file is not included in the main site ZIP download.

None
There is only the prime[] array
primes_db_code.bc
Fixed for February 2013
Code for accessing the packed primes databases below, the main two of which take the form of the standard mathematical nth-prime and prime counting (primepi) functions.
fillprimearray()
Quickly unpacks the first 65,535 primes into the global prime[] array from whichever version of the packed primes database has been loaded at the same time as this file. None of the provided databases contain less than this number so no error should result as long as at least one of them is loaded.
prime(n)
Returns the nth prime; Can only return answers available in whichever version of the packed primes database has been loaded at the same time as this file. Will return an error stating the limitation if it cannot find an answer.
primepi(x)
Returns the number of primes ≤ x; Can only return answers available in whichever version of the packed primes database has been loaded at the same time as this file. Will return an error stating the limitation if it cannot find an answer.
makemods2310_()
Internal function: Creates global metadata arrays used by the above functions for fast access to a packed prime database.
primes_db_minipack.bc
The primes from 13 to 822347 packed using hexadecimal bc code into a global bit array called pd_[].

This and primes_db_code.bc effectively replace the older primes_db.bc database; The old database is well over 1MB in size whereas this and the above are under 53kB.

Like the other data packs below, a global variable, pd_max_ is set so that if data packs are loaded in the wrong order or in a manner such that they overlap, a warning will result. A pd_max_-related load error can mean the database is incomplete or two files have been loaded which contain the same data.

None
There is only the pd_[] array
primes_db_bigpack/0000-0FFF.bc - 1000-1FFF.bc - 2000-2FFF.bc - 3000-3FFF.bc
primes_db_bigpack/4000-4FFF.bc - 5000-5FFF.bc - 6000-6FFF.bc - 7000-7FFF.bc
primes_db_bigpack/8000-8FFF.bc - 9000-9FFF.bc - A000-AFFF.bc - B000-BFFF.bc
primes_db_bigpack/C000-CFFF.bc - D000-DFFF.bc - E000-EFFF.bc - F000-FFFF.bc
The 8.5 million primes from 13 to 151,388,137 packed using hexadecimal bc code into a global bit array called pd_[].

Each consecutive file contains the next 4096 entries of the full database, accessible through the functions in primes_db_code.bc.

In total, these files weigh in at a grand total of 8.53MB, and since bc can take a while to load files given to it on the command line, only as much of the database as is needed can be loaded if the database is split into files this way.

Given that these files take time to load, especially if more than one is used, they display a message after they have loaded so that the user is aware that something is happening rather than a hung bc session.

In each of these packs, a global variable, pd_max_ is set so that if data packs are loaded in the wrong order or in a manner such that they overlap, a warning will result. A pd_max_-related load error can mean the database is incomplete or two files have been loaded which contain the same data.

* These files are not included in the main site ZIP download.

None
There is only the pd_[] array
primes_db_bigpack_onefile/0000-FFFF.bc
The 8.5 million primes from 13 to 151,388,137 packed using hexadecimal bc code into a global bit array called pd_[].

This is the full database available in the above files as one huge block. The data is accessible through the functions in primes_db_code.bc.

This file weighs in at a grand total of 8.53MB, and bc necessarily takes a very long time to load it into memory from disk. During the loading, a message is continually updated on-screen as each block of 4096 entries is loaded into the pd_[] array to avoid giving the impression of a hung bc session. Warning: This file will take on the order of 15 seconds to load on even the fastest 2011 home computer!

In each of these packs, a global variable, pd_max_ is set so that if data packs are loaded in the wrong order or in a manner such that they overlap, a warning will result. A pd_max_-related load error can mean the database is incomplete or two files have been loaded which contain the same data.

* For obvious reasons, this file is not included in the main site ZIP download.

None
There is only the pd_[] array
primes_other.bc
The chaff of interesting but unnecessary functions that would otherwise be included elsewhere. Also, all functions here require funcs.bc in order to run.
aprimepi(x)
An approximation to the prime counting function, π(x); Gives an estimate of the number of primes ≤ x
aq(x)
A Questionable (prime): Turns x into either 2, 3 or base 6 pseudoprime. All primes are of form 2, 3 or 6n±1 for some integer n.
iaq(x)
Inverse of the above; Guaranteed to return a unique value for every prime x.
aq30(x)
A Questionable (prime; base) 30: Turns x into 2, 3, 5 or a base 30 pseudoprime. All primes are of form 2, 3, 5 or 30n±{1,7,11,13} for some integer n.
iaq30(x)
Inverse of the above; Guaranteed to return a unique value for every prime x.
fastguessprime(n)
Uses aprimepi() to find a number very approximately equal to the nth prime. Not guaranteed to return a prime number and almost certainly not the nth prime, but is much faster than the following.
guessprime(n)
Uses the above as well as nearestprime() from primes.bc to find a prime very approximately equal to the nth prime. Almost never actually gives the nth prime, but is guaranteed to return a prime number somewhere near it (usually within 0.5%).
sum_of_factors(x)
Adds together the individual prime factors of a number. e.g. 150 = 2*3*5*5 which becomes 2+3+5+5, and so the result is 15.
sum_of_factor_terms(x)
Adds together the prime power factors of a number. e.g. 150 = 2*3*52 which becomes 2+3+52 = 30.
factor_invert(x)
Raises the powers of prime factors to the power of their primes and re-multiplies
primes_twin.bc
Three functions based upon and requiring the main primes.bc for faster investigation into twin primes.
is_twin_prime(x)
Determines whether x is a member of a prime twin pair. Returns zero if it is not, or else returns the other twin if it is.
nexttwinprime(x)
Returns the nearest twin prime > x
prevtwinprime(x)
Returns the nearest twin prime < x
rand.bc
Modified for February 2013
A floating-point based random number generator. Since bc has no method of obtaining an initial seed, use of randbc is highly recommended. For a sample guessing game which uses this library see guess.bci
rand(x)
For values of x < 0, this will call the srand function with -x as a parameter.
For values of x < 1, will return the previous random number, i.e. the same as the last time the function was called. This is 0 if not previously called.
If x is exactly 1, will return a random floating point number between zero and one (non-inclusive), with as many decimals as the current scale.
For integer values of x > 1, will return an integer between 1 and x.
If x has a floating point portion, this will return an integer between 1 and x+1, with a bias such that floor(x)+1 occurs with the frequency of the fractional part of x.
srand(x)
Resets the internal random seed based on the value of x
thermometer.bc
Temperature scale conversions amongst the five most common scales.
kelvin_to_celcius(k)
kelvin_to_farenheit(k)
kelvin_to_rankine(k)
kelvin_to_reamur(k)
reamur_to_celcius(r)
reamur_to_farenheit(r)
reamur_to_kelvin(r)
reamur_to_rankine(r)
celcius_to_farenheit(c)
celcius_to_kelvin(c)
celcius_to_rankine(c)
celcius_to_reamur(c)
rankine_to_celcius(r)
rankine_to_farenheit(r)
rankine_to_kelvin(r)
rankine_to_reamur(r)
farenheit_to_celcius(f)
farenheit_to_kelvin(f)
farenheit_to_rankine(f)
farenheit_to_reamur(f)
These functions are - hopefully - self explanatory.

General Notes

More information here as and when I think of it...

Disclaimer

Being something of a maths geek, not necessarily on a par with anyone famous or actually intelligent, these files have been created mainly for the fun of it, and also in the hope that someone finds them useful. They have been cobbled together over the years on a sporadic basis. Checks have been made to some degree or another to make sure they work in some fairly unusual conditions, including where scale is set to 0 or where ibase - which affects the interpretation of bc code(!) - has been set to a value other than ten (I would type '10' here, but '10' is the representation of the value of the current number base, whichever base that might be; ten, twelve, sixteen, four hundred and thirty-seven, etc. you see my point.), as well as a reasonable amount of sanity testing.

Before uploading anything to this page, I remove as many bugs as I can find, and hopefully provide sufficient warnings in the above information and the comments in the files themselves about where unusually loose approximations might come about. That said:

I accept no responsibility for any incorrect calculations or otherwise unexpected behaviour, including the consequences thereof, that result(s) from errors that may still exist in these files, whether that be during what would be considered to be sensible use - as part of an interactive bc session, included along with user-written bc code as part of a larger project, etc. - or during any misuse to which they may be put.

Furthermore, I reserve the right to change, update or add to these files at any time, without prior notification, as I see fit.

Contact

If you use these files and find an error, bug, SNAFU, typo or anything else that makes the function you're trying to use behave strangely or just plain wrong, let me know and I'll look into fixing the problem. If you have code suggestions or bug fixes or even questions about bc, feel free to drop me a line.
-- phodd