C Answers to Exercises
Chapter 11
Solution to Exercise 11.1
The normalization ensures that every integer has a unique internal representation, so methods like equals() know that two numbers differ if they don’t have the same number of digits. Without normalization, equals() always has to compare all digits:
public boolean equals(Object other) {
if (!(other instanceof BigNat o))
return false;
for (int i = 0; i < Math.max(size, o.size; i++)) {
if (digitOr0(i) != o.digitOr0(i))
return false;
}
return true;
}
Solution to Exercise 11.2
We need to show that \(c<2\) for the entirety of the algorithm. This is certainly true at the beginning where we initialize \(c\) to 0. Every iteration then updates the carry to \(\lfloor (x_k+y_k+c)/B\rfloor \), which lets us derive an upper bound as follows:
\(\seteqnumber{0}{C.}{7}\)\begin{align*} \bigl \lfloor (x_k+y_k+c)/B\bigr \rfloor &\le \bigl \lfloor \bigl (2(B-1)+c\bigr )/B\bigr \rfloor \\ &< \bigl \lfloor \bigl (2(B-1)+2\bigr )/B\bigr \rfloor &\text {since $c<2$ by induction}\\ &= \lfloor 2B/B\rfloor = 2. \end{align*}
Solution to Exercise 11.3
We need to show that \(-1 \le b \le 0\) for the entirety of the algorithm. This is true at the beginning where we initialize \(b\) to 0.
Under the assumption that the current value of the borrow term satisfies \(-1\le b\le 0\), we derive two bounds for the updated value. For the upper bound we have
\(\seteqnumber{0}{C.}{7}\)\begin{equation*} \bigl \lfloor (x_k-y_k+b)/B\bigr \rfloor \le \bigl \lfloor (B-1 - 0 + 0) / B\bigr \rfloor = 0, \end{equation*}
and for the lower bound
\(\seteqnumber{0}{C.}{7}\)\begin{align*} \bigl \lfloor (x_k-y_k+b)/B\bigr \rfloor &\ge \bigl \lfloor (0 - (B-1) - 1)/B\bigr \rfloor \\ &= \bigl \lfloor -B/B\bigr \rfloor = -1. \end{align*}
Solution to Exercise 11.4
The translation of MultiplyByDigit is shown in Listing C.98.
Listing C.98ch11β―/β―BigNat
public BigNat timesDigit(int y) {
int[] result = new int[size + 1];
long carry = 0;
for (int i = 0; i < size; i++) {
long tmp = (long) result[i]
+ (long) digit(i) * (long) y
+ carry;
result[i] = (int) (tmp & DIGIT_MASK);
carry = tmp >>> DIGIT_BITS;
}
result[size] = (int) carry;
return new BigNat(result);
}
Solution to Exercise 11.5
Let \(c_k\) be the carry at the beginning of the \(k\)th iteration. Since the carry is initialized to 0, we have \(c_0<y\). For every following iteration we have
\(\seteqnumber{0}{C.}{7}\)\begin{align*} c_k &= \lfloor (x_ky+c_{k-1})/B\rfloor \\ &< \lfloor (x_ky + y)/B\rfloor && \text {by induction}\\ &= \left \lfloor \frac {x_k+1}{B}y\right \rfloor \le \left \lfloor \frac {B}{B}y\right \rfloor = y. \end{align*}
Solution to Exercise 11.6
We have \(x<B^n\) and \(y<B^m\), so \(xy<B^nB^m<B^{n+m}\). Any number less than \(B^{n+m}\) can be written using \(n+m\) digits.
Solution to Exercise 11.7
Avogadro’s constant can be computed as follows:
var avogadro = BigNat.fromUnsignedLong(6022140857L);
for (int i = 0; i < 14; i++)
avogadro = avogadro.times(BigNat.TEN);
Alternatively, we can use the pow() method defined in Exercise 11.18:
The digits array of avogadro contains the three base-\(2^{31}\) digits of \(N_A\):
[454574080, 781691483, 130584]
Solution to Exercise 11.8
The value of x is 2. The two calls to slice() first remove the most significant digit \(3\) and then the least significant digit 1, leaving the single digit \(2\).
Solution to Exercise 11.9
No, because in this case the recursion would produce an endless loop.
Solution to Exercise 11.10
In the first iteration we have \(r=0\) and therefore \(t=x_n<B^2\). For each of the following iterations we can use the bound \(r<y\le B-1\) to obtain
\(\seteqnumber{0}{C.}{7}\)\begin{equation*} t = rB + x_i \le (B-1)B + (B-1) = B^2-1 < B^2. \end{equation*}
Solution to Exercise 11.11
For \(x=785, y=18\) we have \(m=2, n=1\), and the algorithm requires two iterations, with a total of 6 correction steps:
\(\seteqnumber{0}{C.}{7}\)\begin{equation*} \begin{array}{ccccc} i & t & \hat q & q_i & r\\ 1 & 078 & \lfloor 7 / 1\rfloor = 7 & 4 & 6\\ 0 & 065 & \lfloor 6 / 1\rfloor = 6 & 3 & 11\\ \end {array} \end{equation*}
Scaling by \(\alpha =\lfloor 9 / (1 + 1)\rfloor =4\) gives us the equivalent problem Divide\((3140, 72)\). We now have \(m=2, n=2\) and the algorithm requires three iterations and no correction steps:
\(\seteqnumber{0}{C.}{7}\)\begin{equation*} \begin{array}{crccc} i & t & \hat q & q_i & r\\ 2 & 031 & \lfloor 03 / 7\rfloor = 0 & 0 & 31\\ 1 & 314 & \lfloor 31 / 7\rfloor = 4 & 4 & 26\\ 0 & 260 & \lfloor 26 / 7\rfloor = 3 & 3 & 44\\ \end {array} \end{equation*}
The true remainder \(44/\alpha =11\) is obtained by dividing the final value of \(r\) by \(\alpha \).
Solution to Exercise 11.12
The implementation of fromDigits() in Listing C.99 is almost identical to the longFromDecimal() method. The main difference is that we are now using a variable radix and that mapping individual digits to their integer values in parseDigit() is slightly more involved.
Listing C.99ch11β―/β―Radix
// Parse integers with a radix between 2 and 36.
public static long fromDigits(String digits, int radix) {
if (radix < 2 || radix > 36)
throw new NumberFormatException("Invalid radix");
if (digits.isEmpty())
throw new NumberFormatException("No digits");
long value = 0;
for (int i = 0; i < digits.length(); i++) {
char c = digits.charAt(i);
value = value * radix + parseDigit(c, radix);
}
return value;
}
private static int parseDigit(char c, int radix) {
int digit = -1;
if (c >= '0' && c <= '9')
digit = c - '0';
else if (c >= 'a' && c <= 'z')
digit = c - 'a' + 10;
else if (c >= 'A' && c <= 'Z')
digit = c - 'A' + 10;
if (digit < 0 || digit >= radix)
throw new NumberFormatException("Invalid digit");
return digit;
}
JDK: NumberFormatException, String
•Radix: parseDigit()
Solution to Exercise 11.13
A possible implementation of toString() is shown in Listing C.100. The outer loop splits the big integer into blocks of 18 digits, and the inner loop converts each block into a string of decimal digits. Since the digits are produced from right to left, the function reverses the list of characters before returning it as a string.
Listing C.100ch11β―/β―BigNat
// Print 'this' as decimal number.
public String toString() {
final int blockSize = 18;
BigNat unit = TEN.pow(blockSize);
StringBuilder sb = new StringBuilder();
BigNat current = this;
do {
var divRem = current.divideAndRemainder(unit);
BigNat next = divRem[0];
long remainder = divRem[1].toUnsignedLong();
for (int i = 0; i < blockSize; i++) {
sb.append(remainder % 10);
remainder /= 10;
if (remainder == 0 && next.equals(ZERO))
break; // omit leading zeros
}
current = next;
} while (!current.equals(ZERO));
return sb.reverse().toString();
}
Solution to Exercise 11.14
The name XYZ27 refers to column \(24\cdot 26^2 + 25\cdot 26 + 26 1 = 16\,900\) and row \(27\). An implementation of CellName is shown in Listing C.101. The fromString() method first computes the column from the sequence of letters at the start of name and then parses the remainder of name as a decimal integer; the toString() method performs the reverse operation. In both cases, the conversion of the column name is similar to a normal radix conversion, except that the digits in this case aren’t zero-based: the numerical value of the smallest digit A is 1.
Listing C.101ch11β―/β―CellName
// Parsing and printing spreadsheet cell names.
public record CellName(int column, int row) {
public static CellName fromString(String name) {
int i = 0;
int column = 0;
while (i < name.length() && name.charAt(i) >= 'A'
&& name.charAt(i) <= 'Z') {
column = column * 26 + (name.charAt(i) - 'A' + 1);
i++;
}
int row = Integer.parseInt(name, i, name.length(), 10);
if (row <= 0)
throw new IllegalArgumentException("Invalid name");
return new CellName(column, row);
}
public String toString() {
var rowStr = new StringBuilder();
int column = this.column;
while (column > 0) {
rowStr.append((char) ('A' + (column - 1) % 26));
column = (column - 1) / 26;
}
return rowStr.reverse().toString() + row;
}
}
JDK: Integer, StringBuilder
Solution to Exercise 11.15
Since computers cannot do arithmetic in either base, we must perform the conversion in two steps. We first convert from base 11 to either a binary (long) or a base-\(2^{31}\) (BigNat) using From-Radix and then convert from this intermediate representation to base 3 using To-Radix.
Solution to Exercise 11.16
If the signs differ, we simply compare them using Integer.compare(); otherwise, we compare their magnitudes and multiply the result by sign (Listing C.102).
Listing C.102ch11β―/β―BigInt
public int compareTo(BigInt y) {
if (sign != y.sign)
return Integer.compare(sign, y.sign);
return sign * magnitude.compareTo(y.magnitude);
}
JDK: Integer
•BigInt
•BigNat: compareTo()
Solution to Exercise 11.17
a)We have \(a=4\) since there are four simpler multiplications and \(b=2\) since each sub-problem involves multiplying numbers with half the number of digits. We therefore have
\[ T(n)=O(n^{\log _2 4})=O(n^2). \]
b)Since only three multiplications are necessary, we have \(a=3\) and \(b=2\). This gives us
\[ T(n)=O(n^{\log _{2} 3})\approx O(n^{1.58}), \]
as stated in the text.
Solution to Exercise 11.18
a)We have \(x^{43}=x^{32}x^{8}x^{2}x\), which gives us three outer multiplications. The individual terms are obtained by repeated squaring
\(\seteqnumber{0}{C.}{7}\)\begin{equation*} x \xrightarrow {\text {square once}} x^2 \xrightarrow {\text {square twice}} x^8 \xrightarrow {\text {square twice}} x^{32}, \end{equation*}
so the total number of multiplications is \(3+5=8\).
b)The pow() function in Listing C.103 computes the integer power of a value of type int by iterating over the binary digits of the exponent.
Listing C.103ch5β―/β―IntMath
// Raise 'value' to the 'n'th power.
public static int pow(int value, int n) {
if (n < 0)
throw new ArithmeticException("Negative exponent");
if (n == 0)
return 1;
int result = 1;
int valuePower = value;
while (n > 1) {
if ((n & 1) == 1)
result *= valuePower;
n >>>= 1;
valuePower *= valuePower;
}
return result * valuePower;
}
JDK: ArithmeticException
c)The number of iterations performed by the while loop in Listing C.103 can be expressed using the bit length of \(k\):
\(\seteqnumber{0}{C.}{7}\)\begin{equation*} \text {\#iterations} = \bitlen (k) - 1 = \lceil \log _2 (k+1)\rceil - 1 \end{equation*}
(See Eq. (9.1) for the definition of the \(\bitlen \) function.) In each iteration, the algorithm performs one multiplication to square the current term and one additional multiplication if the current bit is 1. We can therefore write the total number of multiplications as the sum of the number of iterations, as above, and the number of 1-bits in \(k\), which is the population count of \(k\):
\(\seteqnumber{0}{C.}{7}\)\begin{equation*} \text {\#multiplications} = \lceil \log _2 (k+1)\rceil - 1 + \operatorname {popcount}(k). \end{equation*}
For \(k=43\), this formula evaluates to \(5+4=9\), which is one more than the result obtained in part (a). The additional multiplication is necessary because the first update of result in Listing C.103 is always of the form result = 1 * valuePower.
d)Since the exponent is of type int, we can again use bit operations to iterate over its binary digits. The rest of the algorithm is unchanged, except that all arithmetic operations have to be replaced with the corresponding method calls; see Listing C.104.
Listing C.104ch11β―/β―BigNat
// Raise 'this' to the 'n'th power.
public BigNat pow(int n) {
if (n < 0)
throw new ArithmeticException("Negative exponent");
if (n == 0)
return ONE;
var result = ONE;
var valuePower = this;
while (n > 1) {
if ((n & 1) != 0)
result = result.times(valuePower);
n >>>= 1;
valuePower = valuePower.times(valuePower);
}
return result.times(valuePower);
}
JDK: ArithmeticException
•BigNat: ONE, times()
Solution to Exercise 11.19
a)Every even two-digit number is polydivisible, so there are 45 such numbers; the smallest is 10 and the largest 98.
b)Given a list of all polydivisible numbers with \(k\) digits, all polydivisible numbers with \(k+1\) digits can be computed by appending all possible decimal digits and keeping only those numbers that are evenly divisible by \(k+1\). For an implementation, see Listing C.105.
Listing C.105ch11β―/β―Polydivisible
// Compute a list of all polydivisible numbers.
static List<BigInt> allPolydivisible() {
var result = new ArrayList<BigInt>();
List<BigInt> current = new ArrayList<>();
for (int i = 1; i < 10; i++)
current.add(BigInt.valueOf(i));
for (int k = 1; !current.isEmpty(); k++) {
result.addAll(current);
current = extend(current, BigInt.valueOf(k + 1));
}
return result;
}
static List<BigInt> extend(List<BigInt> current, BigInt newLength) {
var result = new ArrayList<BigInt>();
for (var n : current) {
BigInt nTimes10 = n.times(BigInt.TEN);
for (int digit = 0; digit < 10; digit++) {
var extended = nTimes10.plus(BigInt.valueOf(digit));
if (extended.remainder(newLength).equals(BigInt.ZERO))
result.add(extended);
}
}
return result;
}
There are 20β―456 polydivisible numbers, and the largest one is the 25-digit number:
\(\seteqnumber{0}{C.}{7}\)\begin{equation*} 3\,608\,528\,850\,368\,400\,786\,036\,725. \end{equation*}
The number 3β―816β―547β―290 is the smallest polydivisible number that contains every decimal digit exactly once. Such numbers are also known as pandigital numbers; see [73].
Solution to Exercise 11.20
a)The digits of the integer part of \(n/m\) can be obtained using a simple integer division. We can then compute each digit of the fractional part by computing \(10r/m\), where \(r\) is the remainder of the previous division (Listing C.106).
Listing C.106ch11β―/β―DecimalFraction
// Compute the decimal expansion of x/y to n digits.
public static String decimalFraction(BigInt x, BigInt y, int n) {
var sb = new StringBuilder();
var divRem = x.divideAndRemainder(y);
sb.append(divRem[0].toString()).append('.');
n -= divRem[0].equals(BigInt.ZERO) ? 0 : sb.length() - 1;
for (int i = 0; i < n; i++) {
x = divRem[1];
divRem = x.times(BigInt.TEN).divideAndRemainder(y);
sb.append(divRem[0].toString());
}
return sb.toString();
}
b)The decimal expansion of \(1/998\,001\) starts with all 3-digit numbers between 0 and 999 except for 998 and then repeats itself:
\(\seteqnumber{0}{C.}{7}\)\begin{equation*} 0.000\,001\,002\,003\,004\,005\,006\,007\,008\,009\,010\,\ldots {}997\,999\,000\ldots \end{equation*}
c)As in part (a), we compute the fractional part of the decimal expansion by repeatedly dividing the previous remainder by \(m\). As soon as this remainder takes on one of its previous values, all remainders and digits that follow must repeat themselves. In Listing C.107 we use a hash table remainders to keep track of every remainder that has been encountered so far and the position of the corresponding digit.
Listing C.107ch11β―/β―DecimalFraction
// Compute the decimal expansion of x/y and determine its period.
public static String decimalFraction(BigInt x, BigInt y) {
var divRem = x.divideAndRemainder(y);
if (divRem[1].equals(BigInt.ZERO))
return divRem[0].toString();
var sb = new StringBuilder()
.append(divRem[0].toString()).append('.');
var remainders = new HashMap<BigInt, Integer>();
var rem = divRem[1];
do {
remainders.put(rem, sb.length());
divRem = rem.times(BigInt.TEN).divideAndRemainder(y);
sb.append(divRem[0].toString());
rem = divRem[1];
} while (!remainders.containsKey(rem));
if (!rem.equals(BigInt.ZERO))
sb.insert(remainders.get(rem), "[").append(']');
return sb.toString();
}
The decimal expansion of \(987\,654\,321 / 123\,456\,789\) is
\(\seteqnumber{0}{C.}{7}\)\begin{equation*} 8.[000000072900000663390006036\ldots {}9285243195495713079010989019] \end{equation*}
where the repeating digits are enclosed in square brackets; the number’s period consists of 6β―855β―006 decimal digits.
Solution to Exercise 11.21
a)We use a sign-magnitude representation to store durations: sign holds the sign and digits[] the magnitude in units of seconds, minutes, and so on (Listing C.108). The radix array specifies the radix of each time unit: For seconds and minutes it is \(60\), for hours \(24\), for days \(7\), and for weeks we use an arbitrarily chosen radix of 100β―000. We will use the constants SECOND to WEEK to access to individual array entries.
Listing C.108ch11β―/β―Duration
// A class for computing with time durations.
public class Duration implements Comparable<Duration> {
private final int sign; // -1, 0, +1
private final int[] digits; // units of time
// The radix of each unit of time.
private static final int[] radix = {60, 60, 24, 7, 100_000};
// Constants for different units of time.
private static final int
SECOND = 0, MINUTE = 1, HOUR = 2, DAY = 3, WEEK = 4;
// ...
}
The private constructor in Listing C.109 initializes the sign and digits fields. As in the case of big integers, we normalize durations by requiring that all digits are nonnegative and smaller than the corresponding radix and that the sign is only 0 if all digits are also 0. The implicit assumption is that the digits[] argument points to a newly allocated array that isn’t modified afterward. With this representation, negating durations and computing their absolute value is again trivial.
Listing C.109ch11β―/β―Duration
private Duration(int sign, int[] digits) {
this.digits = digits;
boolean isZero = true;
for (int digit : digits) {
if (digit != 0) {
isZero = false;
break;
}
}
this.sign = isZero ? 0 : sign;
}
public Duration negate() {return new Duration(-sign, digits);}
public Duration abs() {return new Duration(1, digits);}
b)Since durations are stored in normalized form, we can compare (and hash) them by comparing (and hashing) the sign and digits fields (Listing C.110). To implement the compareTo() method, we proceed exactly as for big integers: We first compare the sign of the two durations, and only if they are equal do we compare the digits, from most to least significant.
Listing C.110ch11β―/β―Duration
public boolean equals(Object o) {
return (o instanceof Duration d) && sign == d.sign
&& Arrays.equals(digits, d.digits);
}
public int hashCode() {
return 31 * Integer.hashCode(sign) + Arrays.hashCode(digits);
}
public int compareTo(Duration o) {
if (sign != o.sign)
return Integer.compare(sign, o.sign);
for (int i = radix.length - 1; i >= 0; i--) {
if (digits[i] != o.digits[i])
return sign * Integer.compare(digits[i], o.digits[i]);
}
return 0;
}
c)To construct the duration for a certain unit of time, we set the digits of all shorter units to 0. The remaining digits are computed by taking the remainder v % radix[i] and propagating v / radix[i] to larger units of time (Listing C.111). The listing also shows a method makeSeconds(), which creates a duration for a given number of seconds, and a method seconds(), which returns the current number of seconds. Similar methods for other units of time can be defined in the same way.
Listing C.111ch11β―/β―Duration
// Create a duration measured in a certain unit of time.
private static Duration ofUnit(int value, int unit) {
int sign = Integer.signum(value);
int v = Math.abs(value);
int[] r = new int[radix.length];
for (int i = unit; i < radix.length; i++) {
r[i] = v % radix[i];
v /= radix[i];
}
return new Duration(sign, r);
}
// Create a duration for the given number of seconds.
public static Duration makeSeconds(int s) {return ofUnit(s, SECOND);}
// Return the number of seconds in the current duration.
public int seconds() {return digits[SECOND];}
d)To add and subtract two durations, we use algorithms that are similar to the ones for big integers; the main difference is that each digit now has its own radix. Listing C.112 shows the implementations of addDigits() and subDigits(), which add and subtract the digits of two durations without considering their signs. The subtraction algorithm uses the floorMod() and floorDiv() functions to compute the next digit and the borrow term to correctly handle negative values of tmp.
Listing C.112ch11β―/β―Duration
// Add the digits of two durations, keeping the sign of u.
private static Duration addDigits(Duration u, Duration v) {
int[] result = new int[radix.length];
int carry = 0;
for (int i = 0; i < radix.length; i++) {
int sum = u.digits[i] + v.digits[i] + carry;
result[i] = sum % radix[i];
carry = sum / radix[i];
}
if (carry != 0)
throw new ArithmeticException("Duration out of bounds");
return new Duration(u.sign, result);
}
// Subtract the digits of two durations, keeping the sign of u.
private static Duration subDigits(Duration u, Duration v) {
int[] result = new int[radix.length];
int borrow = 0;
for (int i = 0; i < radix.length; i++) {
int tmp = u.digits[i] - v.digits[i] + borrow;
result[i] = Math.floorMod(tmp, radix[i]);
borrow = Math.floorDiv(tmp, radix[i]);
}
if (borrow != 0)
throw new ArithmeticException("Duration out of bounds");
return new Duration(u.sign, result);
}
JDK: ArithmeticException, Math
•Duration: Duration()
To handle signed durations, we proceed as in the case of signed big integers in Eq. (11.17). The implementations of plus() and minus() are shown in Listing C.113.
Listing C.113ch11β―/β―Duration
// Add two durations.
public Duration plus(Duration v) {
if (sign == v.sign) {
return addDigits(this, v);
} else {
int cmp = abs().compareTo(v.abs());
if (cmp == 0)
return makeSeconds(0);
return cmp > 0 ? subDigits(this, v) : subDigits(v, this);
}
}
// Subtract two durations.
public Duration minus(Duration v) {return plus(v.negate());}
e)To convert a duration to a given unit of time unit, we compute its fractional part (which is determined by digits at positions below unit) and its integer part (which is determined by the remaining digits) and combine the two into a single number of type double. Listing C.114 shows the implementation of totalWeeks(); other methods such as totalDays() are implemented similarly.
Listing C.114ch11β―/β―Duration
private double total(int unit) {
double fracPart = 0.0;
for (int i = 0; i < unit; i++)
fracPart = (fracPart + digits[i]) / radix[i];
double intPart = 0;
for (int i = radix.length - 1; i >= unit; i--)
intPart = (intPart * radix[i]) + digits[i];
return sign * (intPart + fracPart);
}
// Returns the duration as a fractional number of weeks.
public double totalWeeks() {return total(WEEK);}
Solution to Exercise 11.22
a)The implementation of Archimedes algorithm in Listing C.115 starts the iteration with \(p_6\) and \(P_6\). Empirically, this program produces approximately \(0.6\) decimal digits per iteration.
Listing C.115ch11β―/β―PiArchimedes
static double piArchimedes(int iterations) {
double a = Math.sqrt(3) * 2;
double b = 3;
for (int i = 0; i < iterations; i++) {
a = 2 * a * b / (a + b);
b = Math.sqrt(a * b);
}
return (a + b) / 2;
}
Solution to Exercise 11.23
a)
\(\seteqnumber{0}{C.}{7}\)\begin{equation*} \pi \approx 4\left (\frac {4}{5}-\frac {4}{3\cdot 5^3} - \frac {1}{239} + \frac {1}{3\cdot 239^3}\right ) \approx 3.140597 \end{equation*}
b)Machin’s formula tells us that \(100\pi \approx 1600\arctan \frac {1}{5}-400\arctan \frac {1}{239}\). The two arctan functions can be approximated easily since almost all terms are 0:
\(\seteqnumber{0}{C.}{7}\)\begin{align*} 1600\arctan \frac {1}{5} &\approx \left \lfloor \frac {1600}{5}\right \rfloor -\left \lfloor \frac {1600}{375}\right \rfloor =320-4=316\\ 400\arctan \frac {1}{239} &\approx \left \lfloor \frac {400}{239}\right \rfloor = 1 \end{align*} This gives us \(100\pi \approx 316-1=315\).
c)A program that computes approximations to \(U\pi \) is shown in Listing C.116. The computePi() method implements the scaled version of Machin’s formula, and arctan1overX() the scaled version of \(\arctan 1/x\).
Listing C.116ch11β―/β―PiMachin
public static BigInt computePi(BigInt unit) {
BigInt a = unit.times(BigInt.valueOf(16));
BigInt b = unit.times(BigInt.valueOf(4));
return arctan1overX(BigInt.valueOf(5), a)
.minus(arctan1overX(BigInt.valueOf(239), b));
}
static BigInt arctan1overX(BigInt x, BigInt unit) {
BigInt sum = BigInt.ZERO;
BigInt xPower = x;
BigInt xSquared = x.times(x);
int sign = 1;
for (int k = 0; ; k++) {
BigInt denom = BigInt.valueOf(2 * k + 1).times(xPower);
BigInt term = unit.dividedBy(denom);
if (term.equals(BigInt.ZERO))
break;
sum = sum.plus(sign < 0 ? term.negate() : term);
sign *= -1;
xPower = xPower.times(xSquared);
}
return sum;
}
d)Approximating \(10^{10000}\pi \) using Machin’s method produces a 10β―000-digit number that ends in \(5525637566\). The correct digits of \(\pi \) are \(5525637567\), so Machin’s method gives us 9999 correct decimal digits.