Compute pi to base 12 using Python?

[mensanator at aol.com]

Dan Bishop wrote:
> Dick Moores wrote:
> > I need to figure out how to compute pi to base 12, to as many
digits
> as
> > possible. I found this reference,
> <http://mathworld.wolfram.com/Base.html>,
> > but I really don't understand it well enough.
>
> How many stars are in "*************************"?
>
> You probably answered "25".  This means that, for convenience, you've
> broken down the row of stars into ********** + ********** + *****,
that
> is, 2 tens with 5 left over, which the base-10 numeral system denotes
> as "25".
>
> But there's no reason other than tradition why you should arrange
them
> into groups of 10.  For example, you could write it as ******** +
> ******** + ******** + *, or 3 eights plus 1.  In octal (base-8)
> notation, this is written as "31"; the "tens" place in octal
represents
> eights.
>
> In general, in the base-r numeral system, the nth digit to the left
of
> the ones digit represents r^n.  For example, in the binary number
> 11001, the place values for each digit are, right to left, 1, 2, 4,
8,
> and 16, so the number as a whole represents
> 1×16+1×8+0×4+0×2+1×1=16+8+1=25.  This analogous to 25=2×10+5 in
> base-10.
>
> It's also possible to write it as 3×8+0×4+0×2+1×1 = 3001 base 2,
> but by convention, base-r only uses the digits in range(r).  This
> ensures a unique represenation for each number.  This makes "11001"
the
> unique binary representation for decimal 25.
>
> Note that for bases larger than 10, the digits will be numbers that
are
> not single digits in base 10.  By convention, letters are used for
> larger digits: A=10, B=11, C=12, ... Z=35.  For example, the number
> (dec) 2005 = 1×12³+1×12²+11×12+1×1 is represented in base-12 by
> "11B1".
>
> Fractions are handled in a similar manner.  The nth place to the
right
> of the radix point (i.e., the "decimal point", but that term is
> inaccurate for bases other than ten) represents the value
radix**(-n).
>
> For example, in binary,
> 0.1 = 1/2 = dec. 0.5
> 0.01 = 1/4 = dec. 0.25
> 0.11 = 1/2 + 1/4 = 3/4 = dec. 0.75
> 0.001 = 1/8 = dec. 0.125
> 0.01010101... = 1/4 + 1/16 + 1/64 + ... = 1/3
> 0.0001100110011... = 1/10 = dec. 0.1
>
> The last row explains why Python gives:
>
> >>> 0.1
> 0.10000000000000001
>
> Most computers store floating-point numbers in binary, which doesn't
> have a finite representation for one-tenth.  The above result is
> rounded to 53 signficant bits
> (1.100110011001100110011001100110011001100110011010×2^-4), which is
> exactly equivalent to decimal
> 0.1000000000000000055511151231257827021181583404541015625, but gets
> rounded to 17 significant digits for display.
>
> Similarly, in base-12:
>
> 0.1 = 1/12
> 0.14 = 1/12 + 4/144 = 1/9
> 0.16 = 1/12 + 6/144 = 1/8
> 0.2 = 2/12 = 1/6
> 0.3 = 3/12 = 1/4
> 0.4 = 4/12 = 1/3
> 0.6 = 6/12 = 1/2
> 0.8 = 8/12 = 2/3
> 0.9 = 9/12 = 3/4
> 0.A = 10/12 = 5/6
>
> Notice that several common fractions have simpler representations in
> base-12 than in base-10.  For this reason, there are people who
believe
> that base-12 is superior to base-10.
> (http://www.dozenalsociety.org.uk)
>
> > Could someone show me how to do what I need?
>
> You'll need 3 things:
>
> (1) An algorithm for computing approximations of pi.
>
> The simplest one is 4*(1-1/3+1/5-1/7+1/9-1/11+...), which is based on
> the Taylor series expansion of 4*arctan(1).
>
> There are other, faster ways.  Search Google for them.

That one's way too slow. I found the one I use below on Mathworld.

>
> (2) An unlimited-precision numeric representation.  The standard
> "float" isn't good enough: It has only 53 bits of precision, which is
> only enough for 14 base-12 digits.
>
> The "decimal" module will probably work, although of course its
base-10
> internal representation will introduce slight inaccuracies.

I'm using GMPY (see code).

>
> (3) A function for converting numbers to their base-12
representation.

GMPY can do this also.

>
> For integers, this can be done with:
>
> DIGITS = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"
> def itoa(num, radix=10):
>    is_negative = False
>    if num < 0:
>       is_negative = True
>       num = -num
>    digits = []
>    while num >= radix:
>       num, last_digit = divmod(num, radix)
>       digits.append(DIGITS[last_digit])
>    digits.append(DIGITS[num])
>    if is_negative:
>       digits.append("-")
>    digits.reverse()
>    return ''.join(digits)
>
> For a floating-point number x, the representation with d "decimal"
> places count be found by taking the representation of int(round(x *
> radix ** d)) and inserting a "." d places from the right.


# better pi/2 = 1 + 1/3 + (1*2)/(3*5) + (1*2*3)/(3*5*7) + ...

import gmpy
# unlimited precision math library

def pialso(n,b):
# input number of digits (n) in requested base (b)
        p = gmpy.mpq(1,1)
        # gmpy rationals are unlimited precision
        sn = 1
        sd = 3
        c = gmpy.mpq(sn,sd)
        # create the next term to be summed
        num = gmpy.mpf(p.numer())
        # seperately convert the numerator
        den = gmpy.mpf(p.denom())
        # and denominator to gmpy floats
        f = (num/den) * 2
        # to get unlimited precision float
        olds = gmpy.fdigits(f,b,n,0,1,2)
        # extract the requested digits and do base conversion
        done = 0
        while done==0:
                p = p + c
                # sum next term of pi equation
                sn += 1
                # set numerator
                sd += 2
                # and demoniator
                c = c * gmpy.mpq(sn,sd)
                # to create term for next iteration
                num = gmpy.mpf(p.numer())
                # meanwhile, convert this iteration
                den = gmpy.mpf(p.denom())
                f = (num/den) * 2
                # to an unlimited precision float
                s = gmpy.fdigits(f,b,n,0,1,2)
                if s[0]==olds[0]:
                # we're done when the number of digits
                        done = 1
                # we want stops changing
                else:
                # otherwise, keep iterating until we reach the
                        olds = s
                # the desired convergence
        print s[0],sn
        # prints the digits and how many iterations it took
        print

"""

The first 100 digits of pi in base 10.
>>> pialso(100,10)
3141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067
327

The first 100 digits of pi in base 12.
>>> pialso(100,12)
3184809493b918664573a6211bb151551a05729290a7809a492742140a60a55256a0661a03753a3aa54805646880181a3682
353

Note it took more iterations (longer to converge) because base 12
digits are "bigger" than base 10.

"""