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!

If you found this page through an internet search and can't see what you are
looking for, check out the contents box or keep scrolling down.
This page has grown very long and even then some parts are hidden* until they
are activated.

* If you have Javascript disabled, or you're a search engine, everything is visible by
default. This not only allows search engines to classify the page, it allows text-only
HTML renderers like links, lynx, curl and
wget to read the page directly from a terminal - the spiritual home of bc.

Alternatively try the bc FAQ page - you might find
what you are looking for there!

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.

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.

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.

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:

Global variable names and function names which end with underscores indicate that
the object is not intended for use outside the file which contains it.
e.g. output_graph.bc contains a function called or_
which performs a bitwise OR, but is much more simplified than the version found
in logic.bc.

Functions which contain the word "print" in their name write to the console as
well as returning a value. Since bc automatically outputs numbers not assigned
to a variable, it is best to assign these functions to a dummy variable if the
return value is unrequired. See comments in
output_formatting.bc
for an example.

Functions whose names begin with int_ are fast, integer-only
calculations that pay no regard to fractional parts of numbers. Often there
is an equivalent function without the prefix that will work on floating point
numbers. One pair of examples is remainder and
int_remainder which are found in funcs.bc.

Functions whose names begin with is_ determine truth or falsity
and return 1 for true and 0 for false. Many examples of this can be found in
digits.bc.

Functions whose names begin with fast are fast but return
highly approximate results. Often there is an equivalent function lacking the
'fast' in the name which provides better - though not guaranteed perfect -
results, all for the price of a little extra time.

Many functions contain the construct
os=scale; scale=0; /* do something */; scale=os
. This saves the built-in scale variable and switches to
integer-only arithmetic before switching back again. Since bc has no stack to
store system variables, most functions create their own instance of
os, an abbreviation of "old scale" so as not to interfere with
other functions and any interactive session that might be under way.

Fairly often in the code, the constant A appears. This is one of
bc's quirks, in that the capital letters A to F represent the numbers ten to
fifteen regardless of which number base is set within the built-in
ibase variable. Much of bc is designed with base ten in mind,
especially the aforementioned scale variable, which designates
precision in decimal places. Use of A guarantees that
the value obtained is indeed ten and not some other value.

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:

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.

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

Rewritten for June 2012
A suite of functions for basic continued fraction analysis. Uses the
undocumented pass-by-reference feature of arrays where necessary.

Continued fractions

Rational approximation

Upscale / increase rational number accuracy

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.

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 Expansions

Infinite Egyptian Fraction

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).

New for June 2012
Miscellaneous functions that would otherwise be in the main
cf.bc library.

Conway Box

Minkowski Question-mark

Continued fraction transformation

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.

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

Greedy Egyptian fractions

Sylvester Expansions

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).

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()

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.

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.

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.

Digital sum

Reverse

Palindromes

Stringification

Cantor reinterpretation

Bijective (zero-less) numeration

Negative integer bases (negabase)

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.

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.

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.

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.

Counting calculator display segments

Multiply digits

Count digits into arrays

Negapalindromes

Pan digital halving index

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.

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

Combination (Binomial Coefficients)

Permutation

Derangements

LCM factorial / LCMultorial

Euler gamma constant

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.

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.

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

New for February 2013
Near-aliases for some functions in factorial.bc, giving
the functions more classical function names.

Gamma function

Beta function

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))

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.

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.

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.

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)).

Perform numerical integration and differentiation of a single variable function.

Numerical Integration

Numerical Differentiation

Guessing convergence limits

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.

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.

Interest

Savings & Loan

Amortisation / Mortgage

Financial

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.

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).

Fixed word size

Infinite word size

Common bitwise

Twos complement

Bit shifting

Gray code

'Multiplication'

Floating point

Floating point + 'Multiplication'

Gray code + Floating Point

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

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.

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.

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.

Inverse 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.

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.

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.

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.

Bitwise striping

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.

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

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

Prime counting function

Finding the nth prime

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.

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.

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.

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.

The chaff of interesting but unnecessary functions that would otherwise be
included elsewhere. Also, all functions here require funcs.bc
in order to run.

Prime counting function (approximation)

Guessing the nth prime

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

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.

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