- 2016 is the sum of the divisors of both the 42th and the 47th generalized pentagonal number.
- 2016 = 3^3+...+9^3 is the sum of 7 positive cubes.
- The sum of divisors of 2016 squared is a multiple of 2016.
- The 2016th Fibonacci number is a multiple of the sum of divisors of 2016 (= 6552)
- Both 2016 and 2016^2 has 0 as smallest digit and 6 as largest digit.
- The product of the proper divisors of 2016 is a multiple of the sum of the proper divisors of 2016.
- Both the number of divisors of 2016 and the number of positive integers relatively prime to 2016 that are less than or equal to 2016 are perfect squares.
- 2016 = 63*64/2 is a triangular number, i.e, it is the sum 1+2+...+63.
- 2016 = 6!+36^2.
- 2016 (in base 10) is divisible by the sum of their digits in every base from 2 through to 16.
Monday, February 8, 2016
Chinese New Year 2016
Today, February 8, 2016, is Chinese New Year. Just for fun, I looked through OEIS to see what interesting facts I can find about the leap year 2016. Here are some facts about the number 2016:
Monday, December 21, 2015
Another fan awakens
While rewatching Star Wars episodes 1-6, my wife remarked that Han Solo can be roughly translated into Chinese as 單身漢 (Han sounds like 漢 and solo means single which is 單身) which has the meaning of a bachelor, a suitable way to describe his character.
She also wished that she has the Jedi's mind control powers. I told her that she already has that power. For instance, when she says "You will take out the trash" and waves her hand, I take the trash outside to the garbage bin.
She also wished that she has the Jedi's mind control powers. I told her that she already has that power. For instance, when she says "You will take out the trash" and waves her hand, I take the trash outside to the garbage bin.
Thursday, December 17, 2015
The fan awakens
As I mentioned in this post, I have been a longtime fan of Star Wars, as Return of the Jedi was released when I was a teenager. To prepare ourselves to watch the new Star Wars movie, my family and I have been watching the first 6 episodes again over the last couple of weeks. It has been years since I have seen some of them, and it brought back some nice memories. I always chuckled at some of names that are used, for instance Admiral Ackbar who looks like a cuttlefish, is from the species Mon Calamari. Christopher Lee, whose is famous for playing Count Dracula, plays a character called Count Dooku. And midichlorian sounds like mitochondrion, which are microscopic parts in cells responsible for supplying components for generating energy in the cell .
One thing I remember from my childhood is that there is one part of the Star Wars theme song that sounds exactly like part of a deodorant commercial jiggle I saw on a Venezuelan TV Channel.
One thing I remember from my childhood is that there is one part of the Star Wars theme song that sounds exactly like part of a deodorant commercial jiggle I saw on a Venezuelan TV Channel.
Thursday, November 19, 2015
Compression of random data
Over the years, several companies have purported to come up with algorithms that can compress any files, even files of random-looking data. See here for a review of such claims over the years. Compression here means that the compressed file is of a smaller size than the original file. The problem with this is that if you keep applying the algorithm to the compressed file over and over again, you would reach a file of the smallest size: a file of a single bit. Since there are only two possible files of a single bit ("0" or "1"), this simply cannot be the case. This uses the pigeonhole principle and a similar argument of the impossibility of such an algorithm is as follows. There are 2n different bit strings of n bits. On the other hand, there are at most 2n-1 different bit strings of n-1 or less bits. The pigeonhole principle implies it is not possible to compress all n bits strings to a smaller number of bits.
Based on these arguments, Mike Goldman proposed the following challenge (excerpted from comp.compression faq):
This appears to be an easier challenge since the compression and decompression algorithm can be tuned to the data generated. Incidentally, in the following exchange, an attempt is made to win the challenge by exploiting how data is stored in a computer disk system and using the filename as additional storage by splitting the data into many files.
There are several aspects that makes Mike Goldman's challenge possibly solvable and not subject to the constraints of the pigeonhole principle. First of all, the compression algorithm can be tuned to the data. In other words, the compression algorithm only needs to compress only one single file. It doesn't need to compress all files and thus the pigeonhole principle does not apply. Secondly, the challenge did not prescribe what computer language the decompressor can be written in. Can it be written in C, Perl, Python or Matlab? Using a high level language essentially provides the decompressor with complex subroutines that are callable with simple calls, and thus some of the compression can be "hidden" in this way.
Based on these arguments, Mike Goldman proposed the following challenge (excerpted from comp.compression faq):
Mike Goldman <whig@by.net> makes another offer: I will attach a prize of $5,000 to anyone who successfully meets this challenge. First, the contestant will tell me HOW LONG of a data file to generate. Second, I will generate the data file, and send it to the contestant. Last, the contestant will send me a decompressor and a compressed file, which will together total in size less than the original data file, and which will be able to restore the compressed file to the original state. With this offer, you can tune your algorithm to my data. You tell me the parameters of size in advance. All I get to do is arrange the bits within my file according to the dictates of my whim. As a processing fee, I will require an advance deposit of $100 from any contestant. This deposit is 100% refundable if you meet the challenge.
This appears to be an easier challenge since the compression and decompression algorithm can be tuned to the data generated. Incidentally, in the following exchange, an attempt is made to win the challenge by exploiting how data is stored in a computer disk system and using the filename as additional storage by splitting the data into many files.
There are several aspects that makes Mike Goldman's challenge possibly solvable and not subject to the constraints of the pigeonhole principle. First of all, the compression algorithm can be tuned to the data. In other words, the compression algorithm only needs to compress only one single file. It doesn't need to compress all files and thus the pigeonhole principle does not apply. Secondly, the challenge did not prescribe what computer language the decompressor can be written in. Can it be written in C, Perl, Python or Matlab? Using a high level language essentially provides the decompressor with complex subroutines that are callable with simple calls, and thus some of the compression can be "hidden" in this way.
Thursday, October 22, 2015
Back to the future of energy storage
Yesterday was "Back to the Future" day, the day when Michael J. Fox's character arrived into the future from the movie "Back to the Future II". The time machine he used was a DeLorean using a flux capacitor requiring 1.2 Gigawatts of power. I have written a blog post about flux capacitors a couple of years ago and mused that flux capacitors can be considered as normal capacitors and that a large enough supercapacitor today can generate that amount of power if discharged quickly enough. There were endless discussions about whether the prediction about the future in the movie was accurate or not and one comparison to the DeLorean was the Tesla electric car, especially the Tesla Model X with its Falcon wing doors. While the Tesla cannot fly and still needs roads, it does have autopilot capabilities. A fully loaded Tesla comes with a 90kWh battery and if it all discharges within a quarter of a second, it would generate the 1.2 Gigawatts needed for time travel.
Monday, September 28, 2015
Supermoon lunar eclipse
Yesterday was a "supermoon" lunar eclipse, when the full moon is closest to earth and a lunar eclipse is happening at the same time. At the beginning of the evening it was cloudy, but luckily the sky cleared up and we have a great view of the moon outside our garage. We took some nice pictures of the eclipse:
The next supermoon lunar eclipse will occur in 18 years.
The next supermoon lunar eclipse will occur in 18 years.
Tuesday, July 28, 2015
Closest integer square
Given a real number $x\geq 0$, consider the following 2 functions: $f(x)$ is the closest square of an integer to $x$, and $g(x)$ is the integer closest to $\sqrt{x}$.
In other words, $f(x) = y$ where $y = m^2$ for an integer $m$ and $|x-y| \leq |x-n^2|$ for all integers $n$. In case of a tie, i.e. $x$ (resp. $\sqrt{x}$) is midway between two successive squares (resp. integers), we define $f(x)$ (resp. $g(x)$) to be the smallest such square (resp. integer).
The question we like to ask is: is $g(x)^2 = f(x)$ true? For arbitrary real numbers $x$, the answer is no.
The midpoint between two successive integers $a$ and $a+1$ is $a+\frac{1}{2}$ whose square is $a^2+a+\frac{1}{4}$. On the other hand, the midpoint between two successive squares $a^2$ and $(a+1)^2$ is $a^2+a+\frac{1}{2}$. It is interesting to note that the difference between the 2 midpoints is always $\frac{1}{4}$ regardless of what $a$ is (this is true even if $a$ is not an integer).
Graphically this is illustrated as:
This means that for $a^2+a+\frac{1}{4} < x < a^2+a+\frac{1}{2}$, $g(x) = \sqrt{f(x)}+1$. For instance $2.375$ is closer to $1$ than to $4$, but $\sqrt{2.375} = 1.541...$ is closer to $2$ than to $1$.
On the other hand if
$a^2 \leq x < a^2+a+\frac{1}{4}$ or $a^2+a+\frac{1}{2} < x \leq (a+1)^2$,
then $g(x)^2 = f(x)$.
This analysis also shows the following facts:
In other words, $f(x) = y$ where $y = m^2$ for an integer $m$ and $|x-y| \leq |x-n^2|$ for all integers $n$. In case of a tie, i.e. $x$ (resp. $\sqrt{x}$) is midway between two successive squares (resp. integers), we define $f(x)$ (resp. $g(x)$) to be the smallest such square (resp. integer).
The question we like to ask is: is $g(x)^2 = f(x)$ true? For arbitrary real numbers $x$, the answer is no.
The midpoint between two successive integers $a$ and $a+1$ is $a+\frac{1}{2}$ whose square is $a^2+a+\frac{1}{4}$. On the other hand, the midpoint between two successive squares $a^2$ and $(a+1)^2$ is $a^2+a+\frac{1}{2}$. It is interesting to note that the difference between the 2 midpoints is always $\frac{1}{4}$ regardless of what $a$ is (this is true even if $a$ is not an integer).
Graphically this is illustrated as:
This means that for $a^2+a+\frac{1}{4} < x < a^2+a+\frac{1}{2}$, $g(x) = \sqrt{f(x)}+1$. For instance $2.375$ is closer to $1$ than to $4$, but $\sqrt{2.375} = 1.541...$ is closer to $2$ than to $1$.
On the other hand if
$a^2 \leq x < a^2+a+\frac{1}{4}$ or $a^2+a+\frac{1}{2} < x \leq (a+1)^2$,
then $g(x)^2 = f(x)$.
This analysis also shows the following facts:
- If $x$ is an integer, then evaluating $f(x)$ and $g(x)$ do not result in a tie, as the 2 midpoints above are not integers.
- If $x$ is an integer, then $g(x)^2 = f(x)$ as the interval $[a^2+a+\frac{1}{4}, a^2+a+\frac{1}{2}]$ does not contain any integers.
Saturday, January 24, 2015
Trigonometry of regular polgons
In school we learned that the exterior angle of a regular n-sided polygon is $360/n$ degrees or $2\pi/n$ radians. Similarly the interior angle of a regular n-side polygon is $180(n-2)/n$ degrees or $(n-2)\pi/n$ radians. It can be shown that the trigonometry functions (cos, sin, tan, ...) of a rational multiple of $\pi$ is an algebraic number, i.e, a root of a polynomial with integer coefficients. Gauss showed that $n$ being a power of 2 multiplied by distinct Fermat primes (i.e. primes of the form $2^{2^m}+1$) is sufficient to ensure the polygon can be constructed with straightedge and compass and Wantzel showed that this is also necessary (http://mathworld.wolfram.com/TrigonometryAngles.html). This implies that the expression of $\sin(k\pi/n)$, $\cos(k\pi/n)$, etc. can be written explicitly using radicals. Some expressions for these angles are given in (http://en.wikipedia.org/wiki/Exact_trigonometric_constants) for angles where the number of minimal nesting of radicals is 2 or less.
Other values of $n$ for which the minimal nesting of radicals is 2 or less include $n=40$ and $n=120$. For instance for $n=40$, we have:
\[ \sin\left(\frac{\pi}{40}\right) = \left(\frac{1-\sqrt{5}}{8}\right) \sqrt{2+\sqrt{2}} + \frac{\sqrt{2}}{8} \sqrt{2- \sqrt{2}} \sqrt{5+\sqrt{5}} \]
\[ \cos\left(\frac{\pi}{40}\right) = \left(\frac{\sqrt{5}-1}{8}\right) \sqrt{2-\sqrt{2}} + \frac{\sqrt{2}}{8} \sqrt{2+\sqrt{2}} \sqrt{5+\sqrt{5}} \]
\[ \sin\left(\frac{3\pi}{40}\right) = \left(\frac{-\sqrt{5}-1}{8}\right) \sqrt{2-\sqrt{2}} + \frac{\sqrt{2}}{8} \sqrt{2+ \sqrt{2}} \sqrt{5-\sqrt{5}} \]
\[ \cos\left(\frac{3\pi}{40}\right) = \left(\frac{\sqrt{5}+1}{8}\right) \sqrt{2+\sqrt{2}} + \frac{\sqrt{2}}{8} \sqrt{2-\sqrt{2}} \sqrt{5-\sqrt{5}} \]
\[ \sin\left(\frac{7\pi}{40}\right) = \left(\frac{\sqrt{5}+1}{8}\right) \sqrt{2+\sqrt{2}} - \frac{\sqrt{2}}{8} \sqrt{2- \sqrt{2}} \sqrt{5-\sqrt{5}} \]
\[ \cos\left(\frac{7\pi}{40}\right) = \left(\frac{\sqrt{5}+1}{8}\right) \sqrt{2-\sqrt{2}} + \frac{\sqrt{2}}{8} \sqrt{2+\sqrt{2}} \sqrt{5-\sqrt{5}} \]
\[ \sin\left(\frac{9\pi}{40}\right) = \left(\frac{\sqrt{5}-1}{8}\right) \sqrt{2+\sqrt{2}} + \frac{\sqrt{2}}{8} \sqrt{2- \sqrt{2}} \sqrt{5+\sqrt{5}} \]
\[ \cos\left(\frac{9\pi}{40}\right) = \left(\frac{1-\sqrt{5}}{8}\right) \sqrt{2-\sqrt{2}} + \frac{\sqrt{2}}{8} \sqrt{2+\sqrt{2}} \sqrt{5+\sqrt{5}} \]
Other values of $n$ for which the minimal nesting of radicals is 2 or less include $n=40$ and $n=120$. For instance for $n=40$, we have:
\[ \sin\left(\frac{\pi}{40}\right) = \left(\frac{1-\sqrt{5}}{8}\right) \sqrt{2+\sqrt{2}} + \frac{\sqrt{2}}{8} \sqrt{2- \sqrt{2}} \sqrt{5+\sqrt{5}} \]
\[ \cos\left(\frac{\pi}{40}\right) = \left(\frac{\sqrt{5}-1}{8}\right) \sqrt{2-\sqrt{2}} + \frac{\sqrt{2}}{8} \sqrt{2+\sqrt{2}} \sqrt{5+\sqrt{5}} \]
\[ \sin\left(\frac{3\pi}{40}\right) = \left(\frac{-\sqrt{5}-1}{8}\right) \sqrt{2-\sqrt{2}} + \frac{\sqrt{2}}{8} \sqrt{2+ \sqrt{2}} \sqrt{5-\sqrt{5}} \]
\[ \cos\left(\frac{3\pi}{40}\right) = \left(\frac{\sqrt{5}+1}{8}\right) \sqrt{2+\sqrt{2}} + \frac{\sqrt{2}}{8} \sqrt{2-\sqrt{2}} \sqrt{5-\sqrt{5}} \]
\[ \sin\left(\frac{7\pi}{40}\right) = \left(\frac{\sqrt{5}+1}{8}\right) \sqrt{2+\sqrt{2}} - \frac{\sqrt{2}}{8} \sqrt{2- \sqrt{2}} \sqrt{5-\sqrt{5}} \]
\[ \cos\left(\frac{7\pi}{40}\right) = \left(\frac{\sqrt{5}+1}{8}\right) \sqrt{2-\sqrt{2}} + \frac{\sqrt{2}}{8} \sqrt{2+\sqrt{2}} \sqrt{5-\sqrt{5}} \]
\[ \sin\left(\frac{9\pi}{40}\right) = \left(\frac{\sqrt{5}-1}{8}\right) \sqrt{2+\sqrt{2}} + \frac{\sqrt{2}}{8} \sqrt{2- \sqrt{2}} \sqrt{5+\sqrt{5}} \]
\[ \cos\left(\frac{9\pi}{40}\right) = \left(\frac{1-\sqrt{5}}{8}\right) \sqrt{2-\sqrt{2}} + \frac{\sqrt{2}}{8} \sqrt{2+\sqrt{2}} \sqrt{5+\sqrt{5}} \]
\[ \sin\left(\frac{11\pi}{40}\right) = \left(\frac{1-\sqrt{5}}{8}\right) \sqrt{2-\sqrt{2}} + \frac{\sqrt{2}}{8} \sqrt{2+ \sqrt{2}} \sqrt{5+\sqrt{5}} \]
\[ \cos\left(\frac{11\pi}{40}\right) = \left(\frac{\sqrt{5}-1}{8}\right) \sqrt{2+\sqrt{2}} + \frac{\sqrt{2}}{8} \sqrt{2-\sqrt{2}} \sqrt{5+\sqrt{5}} \]
\[ \cos\left(\frac{11\pi}{40}\right) = \left(\frac{\sqrt{5}-1}{8}\right) \sqrt{2+\sqrt{2}} + \frac{\sqrt{2}}{8} \sqrt{2-\sqrt{2}} \sqrt{5+\sqrt{5}} \]
\[ \sin\left(\frac{13\pi}{40}\right) = \left(\frac{1+\sqrt{5}}{8}\right) \sqrt{2-\sqrt{2}} + \frac{\sqrt{2}}{8} \sqrt{2+ \sqrt{2}} \sqrt{5-\sqrt{5}} \]
\[ \cos\left(\frac{13\pi}{40}\right) = \left(\frac{\sqrt{5}+1}{8}\right) \sqrt{2+\sqrt{2}} - \frac{\sqrt{2}}{8} \sqrt{2-\sqrt{2}} \sqrt{5-\sqrt{5}} \]
\[ \sin\left(\frac{17\pi}{40}\right) = \left(\frac{1+\sqrt{5}}{8}\right) \sqrt{2+\sqrt{2}} + \frac{\sqrt{2}}{8} \sqrt{2+ \sqrt{2}} \sqrt{5-\sqrt{5}} \]
\[ \cos\left(\frac{17\pi}{40}\right) = \left(\frac{-\sqrt{5}-1}{8}\right) \sqrt{2-\sqrt{2}} + \frac{\sqrt{2}}{8} \sqrt{2+\sqrt{2}} \sqrt{5-\sqrt{5}} \]
\[ \sin\left(\frac{19\pi}{40}\right) = \left(\frac{\sqrt{5}-1}{8}\right) \sqrt{2-\sqrt{2}} + \frac{\sqrt{2}}{8} \sqrt{2+ \sqrt{2}} \sqrt{5+\sqrt{5}} \]
\[ \cos\left(\frac{19\pi}{40}\right) = \left(\frac{1-\sqrt{5}}{8}\right) \sqrt{2+\sqrt{2}} + \frac{\sqrt{2}}{8} \sqrt{2-\sqrt{2}} \sqrt{5+\sqrt{5}} \]
\[ \cos\left(\frac{13\pi}{40}\right) = \left(\frac{\sqrt{5}+1}{8}\right) \sqrt{2+\sqrt{2}} - \frac{\sqrt{2}}{8} \sqrt{2-\sqrt{2}} \sqrt{5-\sqrt{5}} \]
\[ \sin\left(\frac{17\pi}{40}\right) = \left(\frac{1+\sqrt{5}}{8}\right) \sqrt{2+\sqrt{2}} + \frac{\sqrt{2}}{8} \sqrt{2+ \sqrt{2}} \sqrt{5-\sqrt{5}} \]
\[ \cos\left(\frac{17\pi}{40}\right) = \left(\frac{-\sqrt{5}-1}{8}\right) \sqrt{2-\sqrt{2}} + \frac{\sqrt{2}}{8} \sqrt{2+\sqrt{2}} \sqrt{5-\sqrt{5}} \]
\[ \sin\left(\frac{19\pi}{40}\right) = \left(\frac{\sqrt{5}-1}{8}\right) \sqrt{2-\sqrt{2}} + \frac{\sqrt{2}}{8} \sqrt{2+ \sqrt{2}} \sqrt{5+\sqrt{5}} \]
\[ \cos\left(\frac{19\pi}{40}\right) = \left(\frac{1-\sqrt{5}}{8}\right) \sqrt{2+\sqrt{2}} + \frac{\sqrt{2}}{8} \sqrt{2-\sqrt{2}} \sqrt{5+\sqrt{5}} \]
Monday, November 10, 2014
Evaluating polynomials at regular-space grids
Remember the successive differences approach we learned in school to analyze series of numbers? Consider the set of 4th powers: 0, 1, 16, 81, 256, 625, 1296, 2401, 4096, 6561, ...
Finding the differences between successive terms we get another sequence:
1, 15, 65, 175, 369, 671, 1105, 1695, 2465, ...
Repeating this procedure, we eventually get:
0, 1, 16, 81, 256, 625, 1296, 2401, 4096, 6561, ...
1, 15, 65, 175, 369, 671, 1105, 1695, 2465, ...
14, 50, 110, 194, 302, 434, 590, 770, ...
36, 60, 84, 108, 132, 156, 180, ...
24, 24, 24, 24, 24, 24, ...
Each row is obtained by adding the numbers in the row below to obtain the next term. We notice that the fifth row is constant, i.e. we can obtain the fourth row by simply adding a constant to a previous term. This shows that we can compute each row simply by addition. This approach is valid for any polynomial and is the essence of Nuttall's algorithm (Ref. [1]) to evaluate polynomials on a equally spaced grid. Consider the 3rd order polynomial $p(x) = x^3$ and we want to evaluate it at $x = x_0+n\cdot dx$ for $n=0,1,2, ...$. For simplicity (and without loss of generality, as we can always translate and scale the polynomial), assume that $x_0 = 0$ and $dx = 1$. Using successive differences, we obtain:
$p(n) = p(n-1) + d_2(n)$ where $d_2(n) = 3n^2-3n+1$.
$d_2(n) = d_2(n-1) + d_1(n)$ where $d_1(n) = 6n-6$.
$d_1(n) = d_1(n-1) + d_0(n)$ where $d_0(n) = 6$.
The main idea is that $d_2(n) = p(n)-p(n-1)$ is a polynomial of a lower degree. We can then compute $d_i(n)$ iteratively from $d_i(n-1)$ by adding $d_{i+1}(n)$. Thus $p(n)$ can be computed using the $3$ additions above with the initial conditions $d_0(0) = 6, d_1(0) = -6, d_2(0) = 1, p(0) = 0$.
Given the initial conditions at $n=0$, for a $m$-degree polynomial, this algorithm will require $mk$ additions to compute $p(n)$ for $n = 1, ..., k$.
For the On-Line Encyclopedia of Integer Sequences (OEIS), there are several sequences where such polynomial evaluations are needed. I have written several Python programs to compute these sequences using Nuttall's algorithm. See for example: A161710, A161715, A235810, A138091, A075904, A243298, A007487.
Consider the interpolating polynomial for the divisors of 24 (OEIS A161710):
$\frac{-6n^7 + 154n^6 - 1533n^5 + 7525n^4 - 18879n^3 + 22561n^2 - 7302n+2520}{2520}$
Since the polynomial evaluates to integers for all integers $n$, given the initial conditions, Nuttall's algorithm evaluates this polynomial at integers $n=1,2,...$ using only integer arithmetic, without the need for multiplication or division:
In Ref. [2], it was shown that the initial conditions take some computation and Nuttall's algorithm is more efficient than the direct approach in terms of the number of multiplication operations needed if the number of gridpoints is larger than $o(m^2)$ where $m$ is the degree of the polynomial. However, for the case above the savings in the number of multiplications is more substantial for 2 reasons. First, the initial conditions can be precomputed and hardcoded as in the code above.
Second, in the above example only integer arithmetic is used and the numbers have increasing number of digits as $n\rightarrow\infty$. The initial conditions have smaller number of bits than the numbers that are multiplied in the direct approach to obtain $p(n)$ when $n$ is large and the number of multiplication operations needed on a fixed wordlength computer increases for larger numbers.
In Ref. [1], Nuttall also considered functions of the form $f(x) = e^{p(x)}$ where $p$ is a polynomial. In this case $f(n)/f(n-1) = e^{p(n)-p(n-1)}$, so the successive differences in the polynomial becomes successive ratios in $f$ and the additions become multiplications.
Extensions of this algorithm to polynomials in multiple variables can be found in Ref. [3].
References:
[1] Albert H. Nuttall, "Efficient Evaluation of Polynomials and Exponentials of Polynomails for Equispaced Arguments," IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. ASSP-35, no. 10, pp. 1486-1487, October 1987.
[2] S. C. Datta Roy and Shailey Minocha, "A Note on "Efficient Evaluation of Polynomials and Exponentials of Polynomails for Equispaced Arguments",", IEEE Transactions on Signal Processing, vol. 39, no. 11, pp. 2554-2556, November 1991.
[3] N. K. Bose, "Efficient Evaluation of Multivariate Polynomials Over Multidimensional Grids," IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. ASSP-35, no. 10, pp. 1630-1631, November 1987.
Finding the differences between successive terms we get another sequence:
1, 15, 65, 175, 369, 671, 1105, 1695, 2465, ...
Repeating this procedure, we eventually get:
0, 1, 16, 81, 256, 625, 1296, 2401, 4096, 6561, ...
1, 15, 65, 175, 369, 671, 1105, 1695, 2465, ...
14, 50, 110, 194, 302, 434, 590, 770, ...
36, 60, 84, 108, 132, 156, 180, ...
24, 24, 24, 24, 24, 24, ...
Each row is obtained by adding the numbers in the row below to obtain the next term. We notice that the fifth row is constant, i.e. we can obtain the fourth row by simply adding a constant to a previous term. This shows that we can compute each row simply by addition. This approach is valid for any polynomial and is the essence of Nuttall's algorithm (Ref. [1]) to evaluate polynomials on a equally spaced grid. Consider the 3rd order polynomial $p(x) = x^3$ and we want to evaluate it at $x = x_0+n\cdot dx$ for $n=0,1,2, ...$. For simplicity (and without loss of generality, as we can always translate and scale the polynomial), assume that $x_0 = 0$ and $dx = 1$. Using successive differences, we obtain:
$p(n) = p(n-1) + d_2(n)$ where $d_2(n) = 3n^2-3n+1$.
$d_2(n) = d_2(n-1) + d_1(n)$ where $d_1(n) = 6n-6$.
$d_1(n) = d_1(n-1) + d_0(n)$ where $d_0(n) = 6$.
The main idea is that $d_2(n) = p(n)-p(n-1)$ is a polynomial of a lower degree. We can then compute $d_i(n)$ iteratively from $d_i(n-1)$ by adding $d_{i+1}(n)$. Thus $p(n)$ can be computed using the $3$ additions above with the initial conditions $d_0(0) = 6, d_1(0) = -6, d_2(0) = 1, p(0) = 0$.
Given the initial conditions at $n=0$, for a $m$-degree polynomial, this algorithm will require $mk$ additions to compute $p(n)$ for $n = 1, ..., k$.
For the On-Line Encyclopedia of Integer Sequences (OEIS), there are several sequences where such polynomial evaluations are needed. I have written several Python programs to compute these sequences using Nuttall's algorithm. See for example: A161710, A161715, A235810, A138091, A075904, A243298, A007487.
Consider the interpolating polynomial for the divisors of 24 (OEIS A161710):
$\frac{-6n^7 + 154n^6 - 1533n^5 + 7525n^4 - 18879n^3 + 22561n^2 - 7302n+2520}{2520}$
Since the polynomial evaluates to integers for all integers $n$, given the initial conditions, Nuttall's algorithm evaluates this polynomial at integers $n=1,2,...$ using only integer arithmetic, without the need for multiplication or division:
A161710_list, m = [1], [-12, 80, -223, 333, -281, 127, -23, 1]
for _ in range(1, 10**3):
for i in range(7):
m[i+1] += m[i]
A161710_list.append(m[-1])
|
In Ref. [2], it was shown that the initial conditions take some computation and Nuttall's algorithm is more efficient than the direct approach in terms of the number of multiplication operations needed if the number of gridpoints is larger than $o(m^2)$ where $m$ is the degree of the polynomial. However, for the case above the savings in the number of multiplications is more substantial for 2 reasons. First, the initial conditions can be precomputed and hardcoded as in the code above.
Second, in the above example only integer arithmetic is used and the numbers have increasing number of digits as $n\rightarrow\infty$. The initial conditions have smaller number of bits than the numbers that are multiplied in the direct approach to obtain $p(n)$ when $n$ is large and the number of multiplication operations needed on a fixed wordlength computer increases for larger numbers.
In Ref. [1], Nuttall also considered functions of the form $f(x) = e^{p(x)}$ where $p$ is a polynomial. In this case $f(n)/f(n-1) = e^{p(n)-p(n-1)}$, so the successive differences in the polynomial becomes successive ratios in $f$ and the additions become multiplications.
Extensions of this algorithm to polynomials in multiple variables can be found in Ref. [3].
References:
[1] Albert H. Nuttall, "Efficient Evaluation of Polynomials and Exponentials of Polynomails for Equispaced Arguments," IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. ASSP-35, no. 10, pp. 1486-1487, October 1987.
[2] S. C. Datta Roy and Shailey Minocha, "A Note on "Efficient Evaluation of Polynomials and Exponentials of Polynomails for Equispaced Arguments",", IEEE Transactions on Signal Processing, vol. 39, no. 11, pp. 2554-2556, November 1991.
[3] N. K. Bose, "Efficient Evaluation of Multivariate Polynomials Over Multidimensional Grids," IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. ASSP-35, no. 10, pp. 1630-1631, November 1987.
Thursday, October 2, 2014
IBM Smarter Cities Challenge in Baton Rouge
I am spending 3 weeks in Baton Rouge, Louisiana participating in the IBM Smarter Cities Challenge. You can find the team blog here: http://smartercitieschallenge.wordpress.com/
Here is a view outside of the window of the office that we are working out of, with a view of the Mississippi river and the Horace Wilkinson bridge:
Here is a view outside of the window of the office that we are working out of, with a view of the Mississippi river and the Horace Wilkinson bridge:
Subscribe to:
Posts (Atom)