PARI/GP is an open source computer algebra system I use frequently in my study of binary numbers. It doesn’t manipulate binary numbers directly — input, and most output, is in decimal — so I use it mainly to do the next best thing: calculate with powers of two. Calculations with powers of two are, indirectly, calculations with binary numbers.
PARI/GP is a sophisticated tool, with several components — yet it’s easy to install and use. I use its command shell in particular, the PARI/GP calculator, or gp for short. I will show you how to use simple gp commands to explore binary numbers.
In particular, I’ll show you how to:
- Generate arbitrarily large or small powers of two
- Check if an integer is a positive power of two
- Sum powers of two
- Find the largest power of two contained in a number
- Convert dyadic decimals to dyadic fractions
- Discover properties of binary representations
In the sections that follow, I’ve chosen to show the text log of commands — formatted for display — rather than screenshots of the command shell.
Generate Arbitrarily Large or Small Powers of Two
Nonnegative Powers of Two
gp can generate any arbitrary nonnegative power of two:
? 2^8 %1 = 256 ? 2^64 %2 = 18446744073709551616 ? 2^200 %3 = 1606938044258990275541962092341162602522202993782792835301376
gp can also generate sequences of nonnegative powers of two, from 20 to 2n:
? for(i=0,16,print(2^i)) 1 2 4 8 16 32 64 128 256 512 1024 2048 4096 8192 16384 32768 65536
There are also a few “fun” ways to generate sequences of nonnegative powers of two:
- By computing the divisors of 2n:
? divisors(2^16) %1 = [1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 16384, 32768, 65536]
- By computing the sum of each row, 0 through n, of Pascal’s Triangle:
? for(i=0,16,print(sum(j=0,i,binomial(i,j)))) 1 2 4 8 16 32 64 128 256 512 1024 2048 4096 8192 16384 32768 65536
- By computing the value of Euler’s phi function for 21 through 2n+1:
? for(i=0,16,print(eulerphi(2^(i+1)))) 1 2 4 8 16 32 64 128 256 512 1024 2048 4096 8192 16384 32768 65536
- By computing the order of 5 mod powers of two for 22 through 2n+2:
? for(i=0,16,print(znorder(Mod(5,2^(i+2))))) 1 2 4 8 16 32 64 128 256 512 1024 2048 4096 8192 16384 32768 65536
Negative Powers of Two
gp can generate any arbitrary negative power of two:
- In fraction form:
? 2^-16 %1 = 1/65536 ? 1/2^16 %2 = 1/65536
- In decimal form:
? 2.^-16 %1 = 0.00001525878906250000000000000000 ? 1/2.^16 %2 = 0.00001525878906250000000000000000 ? 0.5^16 %3 = 0.00001525878906250000000000000000
gp can also generate sequences of negative powers of two, from 2-1 to 2-n, in fraction or decimal form:
? for(i=1,8,print(2^-i)) 1/2 1/4 1/8 1/16 1/32 1/64 1/128 1/256 ? for(i=1,8,print(2.^-i)) 0.500000000000000000000000000000 0.250000000000000000000000000000 0.125000000000000000000000000000 0.0625000000000000000000000000000 0.0312500000000000000000000000000 0.0156250000000000000000000000000 0.00781250000000000000000000000000 0.00390625000000000000000000000000
A Note on Precision
When dealing with decimals, you have to make sure you are using the proper precision. The default precision, 28 significant decimal digits, is enough for the examples above — in fact, more than enough, as judged by the trailing 0s.
Here’s an example where the default precision is not enough:
? 2.^-42 %1 = 0.0000000000002273736754432320594787597656
How do I know there’s not enough precision? First of all, the displayed answer has only 40 decimal places, and 2-42 has 42 decimal places (2-n has n decimal places). Also, the displayed answer does not end with the digit ‘5’; all negative powers of two must end in ‘5’.
Here is 2-42 again, printed with 42 significant digits using the ‘\p’ command:
? \p 42 realprecision = 48 significant digits (42 digits displayed) ? 2.^-42 %1 = 0.000000000000227373675443232059478759765625000000000000
(PARI/GP sets its precision based on machine word boundaries, so the number of significant digits used for computation may be greater than the number of significant digits displayed.)
As you can see, we overshot — trailing 0s were added. Even though 2-42 has 42 decimal places, it only has 30 significant digits; that is, 30 digits after the 12 leading 0s.
Here’s how to display 2-42 exactly, without the trailing 0s:
? \p 30 realprecision = 38 significant digits (30 digits displayed) ? 2.^-42 %1 = 0.000000000000227373675443232059478759765625
Check if an Integer is a Positive Power of Two
An integer is a positive power of two if its only prime factor is 2, as demonstrated by these examples:
? 2^64 %1 = 18446744073709551616 ? factor(18446744073709551616) %2 = [2 64] ? 2^64-1 %3 = 18446744073709551615 ? factor(18446744073709551615) %4 = [3 1] [5 1] [17 1] [257 1] [641 1] [65537 1] [6700417 1]
Sum Powers of Two
A binary number represents a sum of powers of two, so you can convert it to decimal simply by performing addition.
Sum Nonnegative Powers of Two
Converting binary integers to decimal is straightforward. For example, 11001002, which represents 26 + 25 + 22, is 10010:
? 2^6 + 2^5 + 2^2 %1 = 100
Sum Negative Powers of Two
Binary fractionals can be finite or infinite; converting finite binary fractionals to decimal is straightforward. For example, 0.1012 is 5/8, or 0.625, in decimal:
? 2^-1 + 2^-3 %1 = 5/8 ? 2.^-1 + 2^-3 %2 = 0.6250000000000000000000000000
0.1111112 is 63/64, or 0.984375, in decimal:
? sum(i=1,6,2^-i) %1 = 63/64 ? sum(i=1,6,2.^-i) %2 = 0.9843750000000000000000000000
Infinite binary fractionals are more involved; they require infinite sums of negative powers of two. To show a simple example, consider 0.012. This represents an infinite geometric series of all negative powers of two except 2-1. It converges to 1/2:
(This is the same summation shown in this site’s header image.)
? suminf(i=2,2^-i) %1 = 0.5000000000000000000000000000
You can see the series converge by using finite sums, adding more terms each time:
? \p 50 ? sum(i=2,10,2.^-i) %1 = 0.49902343750000000000000000000000000000000000000000 ? sum(i=2,20,2.^-i) %2 = 0.49999904632568359375000000000000000000000000000000 ? sum(i=2,30,2.^-i) %3 = 0.49999999906867742538452148437500000000000000000000 ? sum(i=2,40,2.^-i) %4 = 0.49999999999909050529822707176208496093750000000000 ? sum(i=2,50,2.^-i) %5 = 0.49999999999999911182158029987476766109466552734375
Find the Largest Power of Two Contained in a Number
The largest power of two less than or equal to a positive integer n is 2 raised to the integer part of the base 2 logarithm of n; that is, 2floor(log2(n)), or 
. This is used, for example, in the “subtract largest power of two” method of decimal to binary conversion, or in the formula for the solution to the Josephus Problem.
PARI/GP doesn’t compute base 2 logarithms directly, so you have to do a change of base. For example, the largest power of two in 100 is:
? 2^(floor(log(100)/log(2))) %1 = 64
You have to be careful with larger numbers, especially those that are powers of two or close to one. The logarithm function will require greater precision to give the correct result:
? n=2^50 %1 = 1125899906842624 ? 2^(floor(log(n)/log(2))) %2 = 562949953421312 ? \p 100 realprecision = 105 significant digits (100 digits displayed) ? 2^(floor(log(n)/log(2))) %3 = 1125899906842624
(I don’t know if there’s a systematic way to pick the exact precision required — I picked 100 digits arbitrarily, assuming it would be large enough.)
Convert Dyadic Decimals to Dyadic Fractions
A dyadic fraction is a number that can be written in the form a/2n. For our purposes, a < 2n, so it is simply a sum of negative powers of two. For example, 3/8 = 1/4 + 1/8.
A dyadic decimal is what I call a dyadic fraction in decimal form; for example, 0.375. A dyadic decimal can be converted to a dyadic fraction easily: just write it in the form b/10m and reduce to lowest terms. For example:
? 375/1000 %1 = 3/8 ? 375/10^3 %2 = 3/8
PARI/GP reduces fractions to lowest terms automatically, so there’s nothing further to do.
For longer decimals, like those that represent double-precision floating-point numbers, it’s cumbersome — and error prone — to count decimal places by hand. I wrote this gp script to automate the process:
dyadic(decimal_num_string)=
{
decimal_num = eval(decimal_num_string);
decimal_den_exponent = length(decimal_num_string);
decimal_den = 10^decimal_den_exponent;
dyadic_num = numerator(decimal_num/decimal_den);
dyadic_den = denominator(decimal_num/decimal_den);
dyadic_den_exponent = valuation(dyadic_den,2);
print("Decimal: 0.",decimal_num_string);
print("Dyadic fraction: ",dyadic_num,"/2^",dyadic_den_exponent);
}
It takes as argument the string representation of the numerator of the dyadic fraction, which is just the dyadic decimal without the “0.” prefix (a string preserves leading 0s, which must be counted to determine the proper power of 10 denominator!). It computes the denominator that matches the length of the numerator and then divides the two.
The script has the added bonus of displaying the denominator in the form 2n, since large denominators are otherwise not readily identified as powers of two.
Here’s the script in action:
? \r dyadic
? dyadic("375")
Decimal: 0.375
Dyadic fraction: 3/2^3
? 2.^-8
%1 = 0.003906250000000000000000000000
? dyadic("00390625")
Decimal: 0.00390625
Dyadic fraction: 1/2^8
? dyadic("1000000000000000055511151231257827021181583404541015625")
Decimal: 0.1000000000000000055511151231257827021181583404541015625
Dyadic fraction: 3602879701896397/2^55
? dyadic("59999999999999997779553950749686919152736663818359375")
Decimal: 0.59999999999999997779553950749686919152736663818359375
Dyadic fraction: 5404319552844595/2^53
(You might recognize those last two examples as the double-precision floating-point approximations of the decimal values 0.1 and 0.6, respectively.)
Discover Properties of Binary Representations
You can discover properties of the binary representation of a number in two ways: directly, by converting it to binary and viewing the result, or indirectly, by doing calculations related to its underlying binary representation. PARI/GP lets you do both.
Convert Decimal to Binary
PARI/GP can convert integers to binary, although to an awkward format:
? binary(7) %1 = [1, 1, 1] ? binary(2^50-1) %2 = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
PARI/GP can also convert real numbers to binary, to an even more awkward format:
? binary(0.1) %1 = [[0], [0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1]] ? binary(0.6) %2 = [[0], [1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0]]
You can’t tell from the output whether these are finite or infinite expansions, but we know in this case they are infinite: 0.110 is 0.000112, and 0.610 is 0.10012. They are rounded to the nearest dyadic allowed with the given precision.
More bits are printed if you increase the precision with the ‘\p’ command.
Count Number of Trailing Zeros
You can determine the number of trailing zeros in the binary representation of an integer greater than one without first converting it to binary — all you need to do is factor it. If it has no factors of 2, it has no trailing zeros; if it has a factor of 2, it has a number of trailing zeros equal to the exponent to which 2 is raised. For example, 10010 = 11001002 has two trailing zeros:
? factor(100) %1 = [2 2] [5 2]
The valuation function gives the same answer. It prints the exponent of 2 directly:
? valuation(100,2) %1 = 2
Count Number of Bits
The number of bits in the binary representation of an integer n > 0 is 
:
? n=2^32-1 %1 = 4294967295 ? floor(log(n)/log(2))+1 %2 = 32 ? n=2^32 %3 = 4294967296 ? floor(log(n)/log(2))+1 %4 = 33 ? n=100 %5 = 100 ? floor(log(n)/log(2))+1 %6 = 7
(As mentioned above, for large numbers, be sure to set the precision high enough so that the logarithm gives the correct result.)
Count Number of 1 Bits
The number of 1 bits in the binary representation of an integer n > 0 is computed by checking each bit of the binary representation with the bittest() function:
? n=7 %1 = 7 ? sum(i=1,floor(log(n)/log(2))+1,bittest(n,i-1)) %2 = 3 ? n=2^32-1 %3 = 4294967295 ? sum(i=1,floor(log(n)/log(2))+1,bittest(n,i-1)) %4 = 32 ? n=2^32 %5 = 4294967296 ? sum(i=1,floor(log(n)/log(2))+1,bittest(n,i-1)) %6 = 1
(As mentioned above, for large numbers, be sure to set the precision high enough so that the logarithm gives the correct result.)
bittest() works on the underlying binary representation of a number. It numbers bits from right to left, starting at bit 0.
Counting the number of 1 bits is another way to check whether an integer is a power of two — it is a power of two if it has one, and only one, 1 bit.
You can also compute the distribution of 1 bits for the numbers 0 to 2n – 1 by computing row n of Pascal’s Triangle. For example, here are the binary representations of the numbers 0 through 7:
? for(i=0,7,print(binary(i))) [0] [1] [1, 0] [1, 1] [1, 0, 0] [1, 0, 1] [1, 1, 0] [1, 1, 1]
Here is row 3 of Pascal’s Triangle:
? n=3 %1 = 3 ? for(i=0,n,print(binomial(n,i))) 1 3 3 1
This tells us that
- 1 number has 0 1 bits (0)
- 3 numbers have 1 1 bit (1, 10, 100)
- 3 numbers have 2 1 bits (11, 101, 110)
- 1 number has 3 1 bits (111)
Acknowledgements
Thanks to the PARI/GP mailing list, especially Kurt Foster, John Cremona, Bill Allombert, and Karim Belabas, for their ideas and for answering my questions.
I must draw a polar diagram with PARI/GP. I want to write a scrit for it.
Can you help me?
Rosario,
Sorry, I don’t use PARI/GP for much more than I described above. Maybe you can find your answer in the PARI/GP documentation or through the PARI/GP mailing list.
A very useful page – thank you. You mention that binary(0.6)=[[0],[1,0,0,1…]] produces a result with a very awkward format. Putting it mildly! What sort of object is this, please? And can the result be used in further calculations, or turned into a string so that the 1s can be counted, for example?
@Alan,
It is of type “vector”. You can access each “bit” using array notation. Here are some examples I just whipped up:
Convert to string:
Count 1 bits:
(You can also count the number of 1 bits using bittest() as described in the article.)
Thank you! That looks very straightforward. I will try your examples as soon as I have worked out how to save and access scripts in my new PARI – it seems to have changed in the few years since I last used it – either that, or access to the examples folder needs some tweaking…