text:
## CryptoDB

### Craig Gentry

#### Publications

Year
Venue
Title
2022
EUROCRYPT
Non-interactive publicly verifiable secret sharing (PVSS) schemes enables (re-)sharing of secrets in a decentralized setting in the presence of malicious parties. A recently proposed application of PVSS schemes is to enable permissionless proof-of-stake blockchains to keep a secret" via a sequence of committees that share that secret. These committees can use the secret to produce signatures on the blockchain's behalf, or to disclose hidden data conditioned on consensus that some event has occurred. That application needs very large committees with thousands of parties, so the PVSS scheme in use must be efficient enough to support such large committees, in terms of both computation and communication. Yet, previous PVSS schemes have large proofs and/or require many exponentiations over large groups. We present a non-interactive PVSS scheme in which the underlying encryption scheme is based on the learning with errors (LWE) problem. While lattice-based encryption schemes are very fast, they often have long ciphertexts and public keys. We use the following two techniques to conserve bandwidth: First, we adapt the Peikert-Vaikuntanathan-Waters (PVW) encryption scheme to the multi-receiver setting, so that the bulk of the parties' keys is a common random string. The resulting scheme yields $\Omega(1)$ amortized plaintext/ciphertext rate, where concretely the rate is $\approx 1/60$ for 100 parties, $\approx 1/8$ for 1000 parties, and approaching 1/2 as the number of parties grows. Second, we use bulletproofs over a DL-group of order about 256 bits to get compact proofs of correct encryption/decryption of shares. Alternating between the lattice and DL settings is relatively painless, as we equate the LWE modulus with the order of the group. We also show how to reduce the the number of exponentiations in the bulletproofs by applying Johnson-Lindenstrauss-like compression to reduce the dimension of the vectors whose properties must be verified. An implementation of our PVSS with 1000 parties showed that it is feasible even at that size, and should remain so even with one or two order of magnitude increase in the committee size.
2021
EUROCRYPT
This talk is about the past, present and future of fully homomorphic encryption. I will try to convince you that, over the last decade (or so), FHE has made an amazing transformation from proof of concept to a tool that is actually practical and usable. I will review a few homomorphic encryptions schemes (including non-fully homomorphic ones) to demonstrate commonalities in how we start from a homomorphism, and use masking techniques to defeat common attacks (especially linear algebra) on the homomorphism, to obtain semantically secure homomorphic encryption. And, finally, I will try to delineate the search space for fundamentally new FHE schemes that fall outside of the current framework.
2021
CRYPTO
The inherent difficulty of maintaining stateful environments over long periods of time gave rise to the paradigm of serverless computing, where mostly-stateless components are deployed on demand to handle computation tasks, and are teared down once their task is complete. Serverless architecture could offer the added benefit of improved resistance to targeted denial-of-service attacks, by hiding from the attacker the physical machines involved in the protocol until after they complete their work. Realizing such protection, however, requires that the protocol only uses stateless parties, where each party sends only one message and never needs to speaks again. Perhaps the most famous example of this style of protocols is the Nakamoto consensus protocol used in Bitcoin: A peer can win the right to produce the next block by running a local lottery (mining), all while staying covert. Once the right has been won, it is executed by sending a single message. After that, the physical entity never needs to send more messages. We refer to this as the You-Only-Speak-Once (YOSO) property, and initiate the formal study of it within a new model that we call the YOSO model. Our model is centered around the notion of roles, which are stateless parties that can only send a single message. Crucially, our modelling separates the protocol design, that only uses roles, from the role-assignment mechanism, that assigns roles to actual physical entities. This separation enables studying these two aspects separately, and our YOSO model in this work only deals with the protocol-design aspect. We describe several techniques for achieving YOSO MPC; both computational and information theoretic. Our protocols are synchronous and provide guaranteed output delivery (which is important for application domains such as blockchains), assuming honest majority of roles in every time step. We describe a practically efficient computationally-secure protocol, as well as a proof-of-concept information theoretically secure protocol.
2021
TCC
Private information retrieval (PIR) lets a client retrieve an entry from a database without the server learning which entry was retrieved. Here we study a weaker variant that we call random-index PIR (RPIR), where the retrieved index is an output rather than an input of the protocol, and is chosen at random. RPIR is clearly weaker than PIR, but it suffices for some interesting applications and may be realized more efficiently than full-blown PIR. We report here on two lines of work, both tied to RPIR but otherwise largely unrelated. The first line of work studies RPIR as a primitive on its own. Perhaps surprisingly, we show that RPIR is in fact equivalent to PIR when there are no restrictions on the number of communication rounds. On the other hand, RPIR can be implemented in a “noninteractive” setting (with preprocessing), which is clearly impossible for PIR. For two-server RPIR we show a truly noninteractive solution, offering information-theoretic security without any pre-processing. The other line of work, which was the original motivation for our work, uses RPIR to improve on the recent work of Benhamouda et al. (TCC’20) for maintaining secret values on public blockchains. Their solution depends on a method for selecting many random public keys from a PKI while hiding most of the selected keys from an adversary. However, the method they proposed is vulnerable to a double-dipping attack, limiting its resilience. Here we observe that an RPIR protocol, where the client is implemented via secure MPC, can eliminate that vulnerability. We thus get a secrets-on-blockchain protocol (and more generally large-scale MPC), resilient to any fraction f < 1/2 of corrupted parties, resolving the main open problem left from the work of Benhamouda et al. As the client in this solution is implemented via secure MPC, it really brings home the need to make it as efficient as possible. We thus strive to explore whatever efficiency gains we can get by using RPIR rather than PIR. We achieve more gains by using batch RPIR where multiple indexes are retrieved at once. Lastly, we observe that this application can make do with a weaker security guarantee than full RPIR, and show that this weaker variant can be realized even more efficiently. We discuss one protocol in particular, that may be attractive for practical implementations.
2020
TCC
Blockchains are gaining traction and acceptance, not just for cryptocurrencies, but increasingly as an architecture for distributed computing. In this work we seek solutions that allow a \emph{public} blockchain to act as a trusted long-term repository of secret information: Our goal is to deposit a secret with the blockchain, specify how it is to be used (e.g., the conditions under which it is released), and have the blockchain keep the secret and use it only in the specified manner (e.g., release only it once the conditions are met). This simple functionality enables many powerful applications, including signing statements on behalf of the blockchain, using it as the control plane for a storage system, performing decentralized program-obfuscation-as-a-service, and many more. Using proactive secret sharing techniques, we present a scalable solution for implementing this functionality on a public blockchain, in the presence of a mobile adversary controlling a small minority of the participants. The main challenge is that, on the one hand, scalability requires that we use small committees to represent the entire system, but, on the other hand, a mobile adversary may be able to corrupt the entire committee if it is small. For this reason, existing proactive secret sharing solutions are either non-scalable or insecure in our setting. We approach this challenge via "player replaceability", which ensures the committee is anonymous until after it performs its actions. Our main technical contribution is a system that allows sharing and re-sharing of secrets among the members of small dynamic committees, without knowing who they are until after they perform their actions and erase their secrets. Our solution handles a fully mobile adversary corrupting roughly 1/4 of the participants at any time, and is scalable in terms of both the number of parties and the number of time intervals.
2019
TCC
Homomorphic encryption (HE) is often viewed as impractical, both in communication and computation. Here we provide an additively homomorphic encryption scheme based on (ring) LWE with nearly optimal rate ($1-\epsilon$ for any $\epsilon >0$). Moreover, we describe how to compress many Gentry-Sahai-Waters (GSW) ciphertexts (e.g., ciphertexts that may have come from a homomorphic evaluation) into (fewer) high-rate ciphertexts.Using our high-rate HE scheme, we are able for the first time to describe a single-server private information retrieval (PIR) scheme with sufficiently low computational overhead so as to be practical for large databases. Single-server PIR inherently requires the server to perform at least one bit operation per database bit, and we describe a rate-(4/9) scheme with computation which is not so much worse than this inherent lower bound. In fact it is probably less than whole-database AES encryption – specifically about 2.3 mod-q multiplication per database byte, where q is about 50 to 60 bits. Asymptotically, the computational overhead of our PIR scheme is $\tilde{O}(\log \log \mathsf {\lambda }+ \log \log \log N)$, where $\mathsf {\lambda }$ is the security parameter and N is the number of database files, which are assumed to be sufficiently large.
2019
ASIACRYPT
We describe a somewhat homomorphic GSW-like encryption scheme, natively encrypting matrices rather than just single elements. This scheme offers much better performance than existing homomorphic encryption schemes for evaluating encrypted (nondeterministic) finite automata (NFAs). Differently from GSW, we do not know how to reduce the security of this scheme from LWE, instead we reduce it from a stronger assumption, that can be thought of as an inhomogeneous variant of the NTRU assumption. This assumption (that we term iNTRU) may be useful and interesting in its own right, and we examine a few of its properties. We also examine methods to encode regular expressions as NFAs, and in particular explore a new optimization problem, motivated by our application to encrypted NFA evaluation. In this problem, we seek to minimize the number of states in an NFA for a given expression, subject to the constraint on the ambiguity of the NFA.
2018
PKC
We give a framework for trapdoor-permutation-based sequential aggregate signatures (SAS) that unifies and simplifies prior work and leads to new results. The framework is based on ideal ciphers over large domains, which have recently been shown to be realizable in the random oracle model. The basic idea is to replace the random oracle in the full-domain-hash signature scheme with an ideal cipher. Each signer in sequence applies the ideal cipher, keyed by the message, to the output of the previous signer, and then inverts the trapdoor permutation on the result. We obtain different variants of the scheme by varying additional keying material in the ideal cipher and making different assumptions on the trapdoor permutation. In particular, we obtain the first scheme with lazy verification and signature size independent of the number of signers that does not rely on bilinear pairings.Since existing proofs that ideal ciphers over large domains can be realized in the random oracle model are lossy, our schemes do not currently permit practical instantiation parameters at a reasonable security level, and thus we view our contribution as mainly conceptual. However, we are optimistic tighter proofs will be found, at least in our specific application.
2017
EUROCRYPT
2016
TCC
2015
JOFC
2015
TCC
2015
CRYPTO
2014
CRYPTO
2014
CRYPTO
2014
EUROCRYPT
2014
EUROCRYPT
2014
TCC
2013
PKC
2013
TCC
2013
CRYPTO
2013
CRYPTO
2013
ASIACRYPT
2013
EUROCRYPT
2013
EUROCRYPT
2012
EUROCRYPT
2012
CRYPTO
2012
PKC
2011
EUROCRYPT
2010
CRYPTO
2010
CRYPTO
2010
CRYPTO
2010
EUROCRYPT
2010
EUROCRYPT
2009
TCC
2009
EUROCRYPT
2006
CRYPTO
2006
EUROCRYPT
2006
PKC
2005
ASIACRYPT
2005
CRYPTO
2004
ASIACRYPT
2004
CRYPTO
2004
PKC
2003
EUROCRYPT
2003
EUROCRYPT
2002
ASIACRYPT
2002
EUROCRYPT
2001
ASIACRYPT
2001
EUROCRYPT

Crypto 2014
Asiacrypt 2011
Crypto 2011
PKC 2010
Eurocrypt 2009
PKC 2008
Crypto 2008
Asiacrypt 2007
Crypto 2005
CHES 2003
CHES 2002


text:
# Why does e^(pi*i) + 1 = 0

GravitatisVis
It's a really interesting equation. I'm not sure how it works out though. On a similar note, when I raise a whole number to an imaginary non whole number, say 2^2.5i, sine and cosine become involved for some reason. I'm not sure why that happens. Can someone explain this? I finished precalc this year; on to calculus I when school starts in a month. Thanks.

Homework Helper
Gold Member
It's very cute yes. Write e^(pi*i) as cos and i*sin, you'll see why it is true.

Homework Helper
GravitatisVis said:
It's a really interesting equation. I'm not sure how it works out though. On a similar note, when I raise a whole number to an imaginary non whole number, say 2^2.5i, sine and cosine become involved for some reason. I'm not sure why that happens. Can someone explain this? I finished precalc this year; on to calculus I when school starts in a month. Thanks.
So what you need is a definition of e^z=exp(z) that works for complex numbers z. The best motivation for such is analytic continuation, but that requires to much theory, so one of these is often employed.
1) Infinite series
$$exp(z):=\sum_{k=1}^{\infty}\frac{z^k}{k!}$$
2)differential equation
exp'(z)=exp(z) for all complex z
exp(0)=1
3) functional relation
exp(a+b)=exp(a)exp(b) all complex a and b
$$\lim_{z\rightarrow 0}\frac{exp(z)-1}{z}=1$$
4) unmotivated definition
exp(z)=exp(Re(z))cos(Im(z))+i exp(Re(z))sin(Im(z))
5) various ad hoc definitions that state what properties of exp are wanted and define a unique function.

I like 3) best
consider
z=x+i y
exp(z)=exp(x+i y) for x and y real
then
exp(z)=exp(x)exp(i y)
We already know what to do with exp(x) so we turn attention to exp(i y)
exp((x+y)i)=exp(x i)exp(y i)
now for convienance we will define
exp(i x)=C(x)+i S(x)
exp((x+y)i)=[C(x)C(y)-S(x)S(y)]+i[S(x)C(y)+C(x)S(y)]
C(x+y)=C(x)C(y)-S(x)S(y)
S(x+y)=S(x)C(y)+C(x)S(y)
now if x=x y=-x
1=exp(0)=exp(i(x-x))=exp(i x)exp(-i x)
so exp(-i x)=1/exp(i x)=1/(C(x)+i S(x))
=(C(x)-i S(x))/(C(x)^2+S(x)^2)=exp(-i x)/(C(x)^2+S(x)^2)
hence
1=C(x)^2+S(x)^2 for all real x
Now consider
$$\lim_{z\rightarrow 0}\frac{exp(z)-1}{z}=1$$
$$\lim_{z\rightarrow 0}\frac{exp(z)-1}{z}=\lim_{y\rightarrow 0}\frac{C(y)+i S(y)-1}{i y}$$
since we can let z->0 as i y->0
looking at the real part
$$\lim_{y\rightarrow 0}\frac{S(y)}{y}=0$$
Now collecting what we know
C(x+y)=C(x)C(y)-S(x)S(y) for all real x and y
S(x+y)=S(x)C(y)+C(x)S(y) for all real x and y
1=C(x)^2+S(x)^2 for all real x
$$\lim_{y\rightarrow 0}\frac{S(y)}{y}$$
We recognize this can only be if C(x) and S(x) are our beloved functions
C(x)=cos(x)
S(x)=sin(x)
so
exp(x+i y)=exp(x)(cos(x)+i sin(y))
exp(x+i y)=exp(x)cos(x)+i exp(x)sin(y)
for complex exponetials we define
u^v=exp(v log(u))
where the log used is the natural log so log(e)=1
thus the sin and cos in complex exponentials
2^2.5i=exp(2.5 i log(2))
2^2.5i=cos(2.5 log(2))+i sin(2.5 log(2))
2^2.5i~-.161363+.986895i

Last edited:
leon1127
GravitatisVis said:
It's a really interesting equation. I'm not sure how it works out though. On a similar note, when I raise a whole number to an imaginary non whole number, say 2^2.5i, sine and cosine become involved for some reason. I'm not sure why that happens. Can someone explain this? I finished precalc this year; on to calculus I when school starts in a month. Thanks.

son... that is way beyond cal1. It can be proven by series in cal2. But it wont really be decussed until complex analyze.
Proof from MathWorld

Last edited:
Homework Helper
leon1127 said:
son... that is way beyond cal1.
It is a nice precalc/trig type question. The problem is that while many precalc classes cover this material they tend not to be rigorous. Then in introductory calc often due to the large amount of material and the fact that most answers are likely to raise further questions, complex number issues are glossed over. These type of thing tend to be picked up in somewhere along the way and then are delt with in more detail in complex variables or advanced calculus. The idea is simple we want the complex exponential funtion to inherit the sum rule of the real exponential, this along with the rules of complex numbers yeild the famus euler formula.

leon1127
lurflurf said:
It is a nice precalc/trig type question. The problem is that while many precalc classes cover this material they tend not to be rigorous. Then in introductory calc often due to the large amount of material and the fact that most answers are likely to raise further questions, complex number issues are glossed over. These type of thing tend to be picked up in somewhere along the way and then are delt with in more detail in complex variables or advanced calculus. The idea is simple we want the complex exponential funtion to inherit the sum rule of the real exponential, this along with the rules of complex numbers yeild the famus euler formula.

yeah, i agree what you have said. But i dont see it can be proven with pre-cal level math technique. That is why i said it is way beyond pre-cal level. You used limit technique to proof it, but i dont think limit and series is taugh in pre-cal or was taugh in that deep (you were using complex limit and power series, precal barely touch geometric series and real limit).

SteveRives
It may help to back-up a few steps and ask a different question. Namely, what does it mean to raise anything to a complex number? That is, forget for a moment about e and pi related to i. Think instead: What is it to take a thing to the i power? Before we can answer that, we need to figure out the meaning of i. So I ask, What is i in terms of what it does to a graph, or in terms of what it does to a number when it runs into it?

I type with one arm and a baby in hand! So I keep it short and suggestive.

--SR

JonF
Take the McLaurin or Taylor of e^(x*i), i*sin(x) and i*cos(x) and all will be revealed

GravitatisVis
SteveRives said:
It may help to back-up a few steps and ask a different question. Namely, what does it mean to raise anything to a complex number? That is, forget for a moment about e and pi related to i. Think instead: What is it to take a thing to the i power? Before we can answer that, we need to figure out the meaning of i. So I ask, What is i in terms of what it does to a graph, or in terms of what it does to a number when it runs into it?

I type with one arm and a baby in hand! So I keep it short and suggestive.

--SR

That's a good question. I find it very difficult to visualize a number raised to the i power. My knowledge of complex numbers is slim to nil to say the least; i = root(-1). I’m just gonna make a guess here, but does multiplying a number on the complex plane flip it ninety degrees? 4i and 4 form a 90 degrees angle at the origin.

Homework Helper
Gold Member
GravitatisVis:

You asked a simple enough question to answer. I am surprised nobody has given you a simple answer yet.

GravitatisVis said:
On a similar note, when I raise a whole number to an imaginary non whole number, say 2^2.5i, sine and cosine become involved for some reason. I'm not sure why that happens.

It happens because of this definition:

$$e^{ix} \equiv \cos x + i\sin x$$

Using the definition we have:

$$e^{i \pi} = \cos \pi + i \sin \pi = -1 + i \times 0 = -1$$

and therefore

$$e^{i \pi} + 1 = 0$$

If you want something like

$$2^{i \theta}$$

then you need to know that

$$2^x = e^{x \ln 2}$$.

Does that help?

Gza
found this to be funny on the mathworld page:

Gauss is reported to have commented that if this formula was not immediately obvious, the reader would never be a first-class mathematician

Homework Helper
James R said:
GravitatisVis:

You asked a simple enough question to answer. I am surprised nobody has given you a simple answer yet.
Sombody has given a simple answer

It happens because of this definition:

$$e^{ix} \equiv \cos x + i\sin x$$
This is one possible definition. It is a special case of my #4. I call it unmotivated because it is not obvious why any one would think to make that definition. It would also be helpful to know what propertys of the real exponential this preserves and which it does not. Also why is the definition better than another say
exp(i x)=cos(2x)+i sin(2x)

Using the definition we have:

$$e^{i \pi} = \cos \pi + i \sin \pi = -1 + i \times 0 = -1$$
Again the main trouble is of the five commonly used definitions this one gives no reason why it should be used. The other four are more clear, but are difficult to explain on the precalc. level since they use various forms of limits.

Homework Helper
GravitatisVis said:
That's a good question. I find it very difficult to visualize a number raised to the i power. My knowledge of complex numbers is slim to nil to say the least; i = root(-1). I’m just gonna make a guess here, but does multiplying a number on the complex plane flip it ninety degrees? 4i and 4 form a 90 degrees angle at the origin.
Yes this is a way to do it.
consider the unit circle in the complex plain
numbers for which |z|=1
or x^2+y^2=1 (if z=x+i y for x and y real)
the "obvious" way of parametrizing this is
f(t)=cos(t)+i sin(t)
f(s)f(t)=(cos(s)+i sin(s))(cos(t)+i sin(t))
=[cos(s)cos(t)-sin(s)sin(t)]+i[sin(s)cos(t)+cos(s)sin(t)]
=cos(s+t)+i sin(s+t)
=f(s+t)
so
f(s)f(t)=f(s+t)
this looks like an exponential function but it is hard to fit all the pieces together.
f(n t)=f(t)^n
is obvious when n is an integer
with a litle thought it can be seen
f(x t)=f(t)^x
where x is a real number
the trouble is we do not know how to do
f(z t)=f(t)^z
when z is a complex number
this is why since we can define
u^v=exp(v log(u))
for real numbers we take this definition for complex numbers
we wil need to define exp(z) for z complex
for this we can adopt one of the five commonly used (equivelent) definitions of
exp(z) for z complex
at some point we are led to conclude that
exp(a+b)=exp(a)exp(b)
holds when a and b are complex
hence
exp(x+i y)=exp(x)exp(i y)
it turns out we want
exp(i y)=f(y)=cos(y)+i sin(y)
because what we like about exp is
exp(a+b)=exp(a)exp(b)
and
exp(x)~1+x if x is a small number
f(t) is a good canidate for exp(i y) since
f(a+b)=f(a)f(b)
and
f(y)~1+i y if y is a small number

Last edited:
Homework Helper
Gold Member
lurflurf:

Again the main trouble is of the five commonly used definitions this one gives no reason why it should be used. The other four are more clear, but are difficult to explain on the precalc. level

And the initial poster clearly said that he had only just finished precalc.

As a mathematician, it's good to know that you want everything justified and proven from first principles. But somebody who knows nothing about complex numbers is probably better off being presented with some basic results as a first step. They can worry about detailed proofs and motivations later.

Gold Member
I don't suppose there is a non-math model of this formula? Something that demonstrates how these numbers are related in the real (well, or imaginary) world?

SteveRives
GravitatisVis said:
That's a good question. I find it very difficult to visualize a number raised to the i power. My knowledge of complex numbers is slim to nil to say the least; i = root(-1). I’m just gonna make a guess here, but does multiplying a number on the complex plane flip it ninety degrees? 4i and 4 form a 90 degrees angle at the origin.

Let's back up again and think about what it means to take a thing to a power.

Normally, we might think of repeated multiplication when we think of taking something to a power. But there is another way we can think of this. Taking something to a power it a mapping function.

If we only think of taking powers as repeated multiplication, then of course taking to the i power will confuse us! But, if we think of powers as transforming, or mapping, or moving from one set to another, then we might have a better go at it.

Homework Helper
DaveC426913 said:
I don't suppose there is a non-math model of this formula? Something that demonstrates how these numbers are related in the real (well, or imaginary) world?
I can't find anything in a quick google, but there should be something on Joost Burgi either on the net or in a library that provides a geometric interpretation of logarithms. Burgi actually came up with the idea of logarithms before Napier, even though Napier's work was published earlier.

Burgi did work on "prosthaphaeresis" - a method where trigonometry tables were used as an aid in multiplication problems. He expanded on this to develop a geometric approach to logarithms, something that hopefully would be simpler to use for multiplication than trig tables.

Napier also did work in prosthaphaeresis, and developed an algebraic approach to logarithms.

Neither came up with the logarithms used today, but the principles were the same, even if the base was different.

[long and off subject, to be sure, but.....]
One thing I noticed as a side note was that you could develop a slide rule based on prosthaphaeresis that could easily be accurate up to four significant digits if skillfully made. That's kind of interesting since the motivation for developing logarithms was to find a shortcut to long multiplication problems. With a mechanical device to do the calculations, the motivation to progress to logarithms would have been much less. Napier probably wouldn't have invented his "Napier's bones", an invention that extended the ability to perform multiplication to the average business person. Euler probably wouldn't have had Napier's logarithms to read while analyzing infinite series and who knows when the natural log as we know it today would have been invented.

As things turned out, http://www.cee.hw.ac.uk/~greg/calculators/napier/about.html [Broken] revolutionized arithmetic for the average person - except now, instead of getting to actually use sticks for multiplication, elementary school students are expected to memorize them. Brigg's converted Napier's logs to base 10 and called them common logs. The first sophisticated calculating device wound up being based on Brigg's common logs by Edmund Gunter, even though it was kind of a pain to use. William Oughtred turned Gunter's static scales into sliding scales, inventing the first slide rule - a calculating device that was expanded to include trig functions and natural log scales (log log scales) and served as a powerful handful calculator for 350 years until the introduction of the electronic calculator (I doubt the development of natural logs would have been delayed for 350 years, though).

In 1971, you could buy a Post Versalog bamboo slide rule with leather case for about $30 ($146 in today's dollars) - a device that could do multiplication, division, squares, cubes, and square roots, solve quadratic equations, had 8 "log-log" scales to give natural logs for numbers from about .0001 to 20000, could perform vector problems using the trig functions ..... or, you could pay $345 ($1685 in today's dollars) for a Canon pocketronic calculator that could add, subtract, multiply, or divide. In fact, an architect named Paul Hemphill had to finance the first four-function calculator he bought in 1970 - he paid $200 down and paid installments on the remaining$1300 (that's about $7700 in today's dollars). He was still paying off his four-function calculator when Hewlett Packard put out the HP-35 scientific calculator for only$395 - the first hand-held calculator as powerful as a slide rule. Now a days, calculators are so cheap that people can graduate from high school without having to learn anything so complex as 17000*80.

Last edited by a moderator:
GravitatisVis
SteveRives said:
Let's back up again and think about what it means to take a thing to a power.

Normally, we might think of repeated multiplication when we think of taking something to a power. But there is another way we can think of this. Taking something to a power it a mapping function.

If we only think of taking powers as repeated multiplication, then of course taking to the i power will confuse us! But, if we think of powers as transforming, or mapping, or moving from one set to another, then we might have a better go at it.

Mr. Rives, could you explain that in a little more detail? The most we've delved into exponents is, as you said, "repeated multiplication." That's why I find it so difficult to wrap my head around something like a term raised to a non whole number. I've never heard of an exponent "mapping" numbers from one set to another; and to be honest, I have no idea what that even means. Thanks, I appreciate all the answers so far guys.

- Christopher

SteveRives
GravitatisVis said:
Mr. Rives, could you explain that in a little more detail? The most we've delved into exponents is, as you said, "repeated multiplication." That's why I find it so difficult to wrap my head around something like a term raised to a non whole number. I've never heard of an exponent "mapping" numbers from one set to another; and to be honest, I have no idea what that even means. Thanks, I appreciate all the answers so far guys.

- Christopher
Let's start with a simple example. Car squared is not a car that is exponentially bigger. If we apply the idea of squaring to non integers (in this case a car), what do we get? We get gibberish, that's what we get!

Likewise, when we take something to a fractional power, like x to the 1.345 power, we are not talking about a geometric object. Nor are we talking about repeated multiplication. What we are talking about is a function (a mapping) that takes an item and transforms it into a different entity.

Think of cryptography. If you take my English sentence and do a mathematical transform on each letter, the end result will be a sentence that is encoded. The end result is that you map the plain text to code. And that is quite useful. The transformation of each letter is a mapping, and might not correspond to some geometric shape.

Extending the example, taking something to a power is like encrypting plain text, and the log is like the reverse (decrypting the cipher). It is a mapping thing. You are transforming one set to another. It's not something you picture in terms of a graph.

My goal in this short reply is to hint at a different way to think of the exponent idea (so that you are free from thinking only in terms of geometric objects or whole number math). I won't go deeper than this for now, but will wait and see if this made any sense.

Cheers,

Steve Rives

format: plain


text:
# Stattleship! Sport Stats API

I’ve been in contact with the team over at Stattleship. They have a cool API that allows you to get various stats for basketball, football and hockey. I used data from that API to create the following data visualization for their blog. The visualization shows the offensive and special team yards gained by each team remaining in the playoffs. The yardage is totaled for the entire season as well as the one playoff game each team played. I’ve displayed the points off of offensive TDs and special teams scoring, and that score is color coded with to wins and loses. A black background is a win, and a white background is a loss.

# Twitter Sentiment — Penguins VS. Rangers Gm 4

Game 4 of the Penguins-Rangers series featured a brief overtime period that overshadowed the rest of the game as far as tweet volume goes. Rangers fans were more negative at the beginning of the game after the Penguins scored their first goal. Twitter volume picked up for both teams during the overtime period and Rangers fans’ tweets spiked when they won the game and continued throughout the night.

# Twitter Sentiment — Penguins vs. Rangers Gm 3

Unfortunately, my Twitter scraper wasn’t looking for the most viral story of the Penguin’s loss to the Rangers. [I’m not linking to it, but it involves a columnist and the Penguins’ GM.] I was able to get general sentiment over the course of the game. There isn’t too much to analyze. There are more Rangers tweets overall, most likely due to increased interest and a larger market, and I’ve annotated where there was scoring during the game. I’ll probably have a few more updates through out the playoffs.

# 2015 Steelers-Ravens Playoff Twitter Infographics

The Steelers-Ravens playoff game gave me a chance to test out a new analytics server and some of the tools I’ve been working on to make Twitter analysis easy using ad hoc Python scripts. So here goes:

There were a lot of Steelers or Ravens colored emojis, black and gold hearts or buttons and the purple devils. Though for some reasons the ‘crying my eyes out’ emoji is by far the most popular in this collection of tweets. The yellow line represents how many unique tweets there were featuring that emoji. For example, 14 of the same of emoji in one tweet would count for 14 in the blue bar, while it would count for just 1 in the context of the yellow line.

Here’s the hashtag use. The #steelers exceeded the #ravens. This looks cool, but it doesn’t tell you much.

Here’s a bar chart that’s a lot easier to read if you want the information.

# Getting Lucky in a Playoff Series

Sports have a constant uncertainty and randomness in every aspect of the game including determining champions. This is one area you wouldn’t expect to have a lot of variability, since you would want the team that has the best roster composition and played the hardest to win the championship. This concept is usually brought up in the arguments against the one-game Wild Card round that MLB introduced in 2012 saying there’s too much that can happen in one game to determine the fate of a season. [The counter-argument to this specifically is the division winners now have a reward for winning the division, besides having cool sweatshirts.]

The Sports Side

The basis for championship series in MLB, NBA, and NHL is the an odd-numbered of games series with the champion being the team that wins the majority of games in those series. Most sports use a 7-game series; so for example the Boston Red Sox had to win 4 games to win the World Series last year. Using an example of randomness I got from Leonard Mlodinow’s The Drunkard’s Walk: How Randomness Rules Our Lives, I can illustrate how a team that’s clearly an underdog can win a playoff series against a superior opponent. Mlodinow has a recorded lecture where he explains what he wrote in his book. [It’s a good book, you should read it.]

Let’s use two teams; one is the Favorite, and it is assumed they will beat the Underdog 55% of the time [given enough games]. This also means that the Underdog will win 45% of the time. This represents win probabilities more uneven that you are likely to find in a playoff game since teams are typically much more evenly matched [at least in baseball]. The last assumption of this example is that the teams win probabilities don’t change with a different starting pitcher or home field/court/ice advantage. These are terrible assumptions if you wanted to project real playoff series, but the underlying principle of random sequencing will still hold true.

In order to win the playoff series, a team has to win a certain number of games before the opponent wins that number. To model this distribution based on pure randomness, you can use the negative binomial distribution to determine the probability that the Favorite will win a 7-game playoff series in 4, 5, 6, or 7 games. If you wanted to design a playoff series to minimize the chances that the underdog will win, you’d want to choose a number of games which would have the smallest probability of the underdog winning the series.

This chart shows the probability for all 8 possible outcomes of a 7-game playoff series based of a 55/45 winning percentage split and pure randomness with no home-field advantage. As you can see there’s a substantial chance [39% probability] that the Underdog will win a 7-game series. 39% is rather large, and this is a 7-game series. Baseball also employs a 5-game series for their division series (LDS) and a one-game playoff for the the Wild Card round (WC). The chances of an upset becomes more likely as the number of games decreases. I’ve also added another set of teams (60/40 split — greater disparity) for comparison’s sake.

It should be obvious that the 1-game series has the greatest chance of an upset, hence the objections to its use in baseball. Though my contention would be that a 3-game series does not offer much more certainty that the best team will win.

The Math Side

I first calculated these probabilities by writing out all the possible combinations then adding up those probabilities. I have since realized there was a much easier way to determine these probabilities, and that is by using the negative binomial distribution (NBD). If you want to familiarize yourself with what the distribution represents please read the count data primer. In short, the NBD will determine the probability that a team will lose a certain number of games [0-3] before the other team wins 4 games. The NBD is defined by the following function:

$P(X=k) = {{r+k-1}\choose{k}} p^{r} (1-p)^{k}$

where X is the random variable whose probability we are calculating, k is number of Team A losses [this will vary], r is the number of Team A wins [for the 7-game series, it will be 4 games], and p is the probability of Team A winning. In the case of this example we will be determining the probability of Team A winning a 7-game series, when Team A has a 55%/45% advantage over Team B.

$P(X=2) = {{4+2-1}\choose{2}} 0.55^{4} (1-0.55)^{2}$

$P(X=2) = {{4+2-1}\choose{2}} 0.55^{4} (1-0.55)^{2}$

$P(X=2) = 10 * 0.0915 * 0.2025 = 18.53\%$

This is the probability for just one possible outcome, Team A wins the series in 6 games. To determine the probability that Team A wins the series, you add the probabilities for each outcome Team A wins in 4, 5, 6, or 7 games. So this calculation then repeated for every loss possibility:

$P(WinningSeries) = P(X=0) + P(X=1) + P(X=2) + P(X=3)$

$P(WinningSeries) = 9.15\% + 16.47\% + 18.53\% + 16.68\% = 60.83\%$

From these calculations, there is a 60.83% chance that the Team A wins just by randomness. Conversely, there is a 39.17% [100% – 60.83%] chance that Team B, the inferior team, wins because of random sequencing.

Conclusion

The MLB Wild Card game rightfully gets criticized for being too susceptible to having a bad day or getting a bad bounce. I wanted to illustrate that any playoff series has a lot of randomness in it. Beyond the numbers, people remember the bad bounces way more than they remember the positive or neutral events that occur [negativity bias]. A bad bounce or a pitcher having a bad day could easily benefit the team you are rooting for. The only real way to root out the randomness you would need to play hundreds of games, and somehow I don’t think that is feasible.

# Twitter Analysis – Penguins Game 7

I’ve been listening to 93.7 The Fan while running the analysis for this, and I never realized that people can say the same thing over and over again but in slightly different ways. Also all tweets were captured AFTER THE CONCLUSION OF THE 1st PERIOD.

Everyone knows Twitter is the best venue to vent your anger about sports teams. I was able to the statistical programming language R to scrape tweets which had certain keywords or hashtags in them, put them in a database and then flag the tweets that contain certain keywords or collection of words. I had about 20 keywords including: “penguins”, “pens”, “rangers”, “game 7”, and “firebylsma”. I also search for any of the handles of the local hockey writers, because a lot of people will reply to the sports writers during the game with their own opinions.

The first graph has the total number of tweets that I scraped and tweets that I flagged as ‘swearing’. For the most part, I feel like if someone swore in the tweet, it indicates anger or at the very least aggressiveness. As mentioned before these graphs begin after the 1st period ends. And tweets containing the keyword ‘rangers’ has been filtered from this first graph to include a greater amount of Penguins fan.

The quickest and most basic analysis is the number of tweets as a time-series. As soon as you look at the time-series line graph, you can tell when the game ends [9:41 PM]. It’s like Mt. Everest in the graph. Looking closer you can see the spikes where each team scores. An interesting occurrence happened right before the end of the game. Twitter got quiet. The tweets per minute dropped below 500 right before it exploded to a few thousand per minute. I am attributing this silence to people actually watching the game during the tense last minute.

The tweets peaked about 2 minutes after the game ended indicating a minimal lag which includes the time of picking up the smart phone, unlocking it, and composing the tweet. However, the angry, swearing tweets peaked right as the game ended indicating more visceral emotions instead of more thought-out, 140-character commentary. There are two severe dips that I can’t account for at 9:49 PM and 10:04 PM. If anyone knows something that occurred at this time, please let me know. Since there is a clear downward trend and the game was no longer being played, I am going to write those dips off as some technical difficulties that didn’t allow a lot of tweets to be sent at those times.

Since Twitter isn’t an invention specific to just Penguins fans, I separated and compared two sets of tweets: one containing the word “penguins” and one comparing “rangers”. There are a lot more Rangers fans than Penguins fans, because the “rangers” tweets outnumber “penguins” tweets at almost every time. The “penguins” tweets did spike when they scored their only goal. Interestingly enough, there was a lot of swearing right when the goal was scored. Penguins fans are just so angry!

Bottom-line calm down; it’s just sports.

# Against All Odds — Upsets

It’s a great time to be a Dayton fan!  It’s the first time the school has reached the Sweet Sixteen since before all the Dayton fans I know where born…and they did it as an 11 seed!  Their game against Ohio State was the first game of the tournament to tip.  A little over two hours later, everyone’s brackets were busted.  After watching the tournament this weekend, I felt like there were a lot of upsets this year.   (Or at the very least my bracket was getting busted up pretty quickly.)  But are there really more upsets this year than normal?

First, I’m going to define an upset as any lower seed beating a higher seed.  I’m of the personal belief that 8/9, 4/5, and 1/2 match-ups shouldn’t count as upsets, but for this analysis, I’m going to consider these as possible upsets.  First, let’s look at how many upsets there were this year.  Through two rounds, there have been 13 upsets.  That’s one less than last year at this time, and just at the average (if you round).  So this is a rather average year.  Three of four 1-seeds are still alive — not too much different from what you might expect.

NCAA Tournament Upsets By Year

Historically, 1999 had the most upsets with 19 in the first weekend of play.  Nothing really stuck out like how Florida Gulf Coast got to the Sweet Sixteen as a 15-seed last year.  1991 had the fewest number of upsets in the first weekend with just nine.  All the upsets through out the years appear to be random noise fluctuating around an average of 12.8 upsets (out of 48 games played) in the first two rounds per year.  The conclusion you can draw from this is that the number of upsets is rather consistent over the years with not much systematic change from year to year.

Thirteen upsets is a lot; it’s almost 1/3 of all the games played this weekend.  Last week, I posted the probability that a seed would win in the first round of the tournament.  This was a linear relationship starting with an almost certain probability for the 1-seeds and then going to a 50/50 split for an 8/9 game.  On the surface it doesn’t seem like it almost 1/3 of the games would be upsets, but if you look at all the possibilities it will make more sense.

Let’s look at Dayton’s 11-seed.  A 11-seed has a historical 34% chance of upsetting a 6-seed in the first round, but when considering there are four distinct 6-11 seed match-ups each year there’s only a 19% chance that all 6-seeds will win their first round games.  In fact, the most likely scenario is that just one 6-seed will upset a 11 seed.  This year there were two 6-11 upsets which is the second most likely scenario at 30% (still more likely than not getting any upsets).

The following table depicts the probability of different scenarios for each first round seeding combination.  All the green area on the table is why everyone’s brackets bust every year.  Keep reading if you are interested in the math, otherwise you might want to bounce, because it’s gonna get boring.

NCAA Tournament First Round Win Probability by Seed

Still here?  Ok.  The basis for determining the probability of the upset scenario is the binomial distribution.  A binomial distribution requires two things, a binary outcome (hence the bi- prefix)  and a set probability of how that outcome is achieved.  The simplest example of a binomial distribution is determining the probability of successive coin flips.  The probability function is given as

$P(X) = (^n_k) p^k q^{1-k}$

The $(^n_k)$ term is the combination of n terms taken k-at a time.
$p$ is the probability of the event happening — the win probability
$q$ is the compliment of the event so in this cause it would be probability of losing
$n$ will be 4 since there are four games for a seed match up
$k$ will be 0-4 depending on how many upsets we are looking for.

Looking at the probability that two (and only two) 11-seeds upset 6-seeds that will be

$P(X) = (^4_2) (.34)^2 (1-.34)^2 = 6 * .34^2 * .66^2 = .302 = 30$%

You can derive this equation by writing out probability trees (if you remember those from high school math).  The problem with that method is that for each outcome (# of upsets = [0, 1, 2, 3, 4]) you have to write out the different combinations of games for each outcome.   This can get unwieldy quickly.  Binomial distributions can be used for many different applications, including the aforementioned coin-flip, likelihood of combinations of boy/girl babies, the probability that the ‘better’ team loses a 7-game playoff series, the likely number of winners for the lottery…so this will rear it’s head again for NHL, NBA, or MLB playoffs.

# 2014 NCAA Tournament Predictions — Monte Carlo

The process of simulating the NCAA tournaments involves two-steps.  The first is determining what statistical prediction model to use to determine the outcome of a game.  The second step is to simulate the entire tournament.   Simulating the tournament multiple times and keeping track of each outcome is called a Monte Carlo simulation. This simulates the entire tournament 10,000 times and tabulates the results from each round.

On to what the computer says [entire bracket png]!  Surprise, surprise the computer says almost everything that you might surmise by using your gut.  It predicts Florida winning  the entire tournament with almost a 20% probability, and it also predicts very few upsets.  Everything is pretty much what you would expect.  As I said earlier in the week, the committee does a pretty good job seeding everybody overall.

There are a few ‘undervalued’ teams the simulation has picked.  Villanova, a 2-seed , is projected to go to the Final Four.  All the other Final Four teams are 1-seeds, which is the seed with the highest probability to reach the Final Four historically.  The most dramatic prediction, I think, is the North Dakota State upset.  NDSt is a 12-seed.  12-seeds are the most frequently undervalued seed performing well above the expected winning percentage.  In betting on an upset, I’d be looking for a lower seed’s win probability to be higher than average for that seed.  However, NDSt’s win probability is not only higher than average, it’s higher than Oklahoma’s win probability.  I’m putting NDSt to win the first game in all my brackets today.

The rest is pretty obvious, but I’m going to be interested to see how this bracket does against what really happens.   And how simulations look after each round.

The bracket is broken into its four regions for viewing ease.  The winners are in green and the upset winners are in yellow.  The Final Four is the last graphic.

# NCAA Tournament — Seeding

All of the analysis is only looking at the 64 team field from the 1985 tournament through 2013.  Before 1985 there were less than 64 teams invited.  Opening Round games are also ignored.

The NCAA Selection Committee just released the seeding for the upcoming tournament, and everyone over the next few days will be filling out brackets to win 50\$ from their office pool.  But who to pick?  Later this week, I’m going to look at a basic way to predict the point spread for an NCAA tournament game.  Before I look at that prediction, let’s examine how the different seeds perform in the tournament.

Seeds are ordinal ranks (1 through 16) placed on the 64 teams and then those teams are pitted against each so the highest seeds won’t face each other till the final round in the Regional Championship (Elite Eight).  So how well do the seeds perform?  The graphs above show how well each seed has done since 1985, the first year of the 64 team bracket.  The 1 seeds are clearly are the most successful in any around, and beyond that there’s a linear pattern with first round wins.  So the better the seed the more likely the team is to go further.

In fact the likelihood of the highest seed winning the first game goes from 100% for the 1-16 game to a coin-flip for the 8-9 game.

Looking at the first round wins, there are two interesting seeds, the 9 and 12 seeds, which outperform the pattern.   This most likely happens because the Selection Committee undervalues teams it assigns to those slots.  There’s an asymmetry of information available for power conferences vs mid-majors.  Mid-majors have less games against top 50 teams, so the the committee has less information to judge the team.

Now as far as selecting teams to win tournament games.  The safest pick is to always pick the highest seed.  The argument against this is that its very boring and you might not win.  Being the most likely to win does not really translate into winning in singular instances.  But the data points to 1 seeds being the overwhelming favorite to get into the Final Four over any other particular seed by a 2:1 margin.

Picking upsets.  The underlying principle of the upset in this tournament is the small sample size of just one game.  A lot can happen in one game.  Even if a team has a 78% chance of winning, which is rather favorable (Like a 4 seed)….there’s still a 22% chance that team will lose.  With 64 teams and 63 games there’s a lot of room for upsets to happen by random chance.  The randomness is in how the team feels that day, how fouls are called, if a foul was called in a critical junction, or come down to how a ball bounces or falls.

Given the way the seeds have performed since the tournament expanded in 1985 to a 64-team tournament.  The selection committee does a good job seeding the teams, since they are following a pattern you would except from seeding the best teams against the worst.

# 2013 AFC Playoffs and Bayesian Statistics

How does a team like the Steelers go from 0-4 to a dark horse for the final AFC playoff spot to inches away from clinching to the birth?  Then the Chargers who were a dark horse themselves go on an secure the 6th seed?

Philip Rivers mentioned, in an endorphin-high interview, that no one gave them a chance, saying the odds were against the Chargers.  This delves into the different realms of motivational and analytical that I want to stay away from, but how did these odds really work and how did they change through out the day?  Week 17 of the 2013 NFL season was a great example of Bayes’ theorem.

Here’s the math for Bayes’ theorem:

This is read: the probability of A given that B happens is equal to the probability of B given that A happens times the probability of A divided by the probability of B.  We are going to use this basic form of Bayes’ theorem to understand what happened to the Steelers and Chargers playoff odds throughout the day on December 29th.

The important point I’m going to illustrate here is that probability is not a property inherent in an event.  It is rather a guess or calculation based on known facts or frame of reference at the time.   When that frame of reference changes (in this case when the NFL schedule plays out) we can calculate new odds.  I’ll be using the above mathematical formula to calculate new probabilities as we discover new information through out the day for the Steelers.  Then I’ll compare this with a similar chart for the Chargers, who made the playoffs.

The first and largest problem for this exercise is determining the win/loss probability for the games.  I’ve searched online and there is much disagreement on what the playoffs odds where at the beginning of the day, let alone what any singular games odds were.  I could use Vegas odds as crowd-sourced odds, however, I’ll make this part easy and just make up numbers for illustrative purposes.  I have the Steelers and Chargers win probability weighted high because they were playing  the Browns and the Chiefs’ back-ups.  The Ravens and Dolphins I have at even odds of 0.5.  These could be endlessly debated, but let’s just assume they are correct.

Assumed probability chart for week 17.

For the Steelers to make the playoffs four things needed to happen.  First they had to win. [P(SteelersW) = .85]  Next the Ravens and Dolphins had to lose . [P(DolphinsL) = P(RavensL) = .5]  And then San Diego had to lose later in the day.  [P(ChargersL) = .75]  If you multiply all these together you will get roughly P(SteelersPlayoffs) = .05.  So there’s a 5% chance with the information at the beginning of the day that the Steelers will make the playoffs.

So let’s use Bayes’ theorem to calculate what the playoff are are if we know the Dolphins lose their game — P(SteelersPlayoffs | DolphinsL).

From just knowing the Dolphins losing their game, you can infer that the Steelers chances of gaining of a playoff birth is twice as likely as it was before.  I should also explain why P(DolphinsLoss | SteelersPlayoffs) is equal to 1.   This term assumes the Steelers made the playoff and asks what the probability is of the Dolphins’ loss given this information.  The Dolphins must lose if the Steelers make the playoffs so the term is equal to 1.  This will be true for every game we are considering.  (This problem becomes more complicated if there are multiple paths to playoffs, because the term will no longer be 1.)

Starting at P(SteelersPlayoffs) = .05, you can calculate the conditional probability in a chain as the Steelers’ win, Dolphins’ loss, and Ravens’ loss occurs.  Then the probability of the Steelers making the playoffs is calculated solely on the remaining game Chiefs/Chargers.  In-game win probabilities are calculated on advancednflstats.com.  They are dependent on time left in the game, score, and field position, and down, and they are independent of team skill.  This is the probability graph for the Chargers game.  I used this for the final two Steelers calculations: just before the Chiefs missed a FG, and then in OT during the Chief’s final drive after the Chargers kicked a FG to take the lead.

Steelers Playoff Probability — Area chart

You can see how the Steelers’ probability changed after each event, and how small the area was until the Chargers almost lost their game.  I will compare this with the Chargers, who had a better chance all through the day, except when the Chiefs were threatening to win.  The large green area is the probability at the end of the Chargers game when they won, and it’s value is 1.0 denoting they have clinched the playoff birth.

Chargers Playoff Probability — Area chart

To respond to Philip River’s on field comments about the Chargers having long odds, it’s misleading to think the Chargers somehow overcame those odds themselves, when they were the favorites to capture the final wild card spot when their game started.  It was the others teams’ losses that increased  the odds in the Chargers favor before they even played.

This exercise serves as an example of how and why probabilities change over time, and it illustrates how probability relies on known information or a reference point.  And how changing the reference point affects the known probability.


text:
# Patterns

My kid’s math teacher recently posted a puzzle for the parents that got me down a rabbit hole of sorts and made me reflect on a number of things. This post is about that … and I mix the math and programming bits. I describe the puzzle in the first section and then the rest reveal aspects of the solution. So if you want to avoid the spoilers, you should stop reading after the first section.

## The Puzzle

Consider an $$N \times N$$ square grid, where each cell may contain a letter-digit combination - like $$A1, A2, B4, D5$$ and so on. The letters are to be drawn from $$ABCDEFGHI$$ and the digits from $$123456789$$. For the purpose of this discussion, we’ll limit ourselves to $$N \leq 9$$, but you can consider larger grids on your own by extension.

The task is to fill the grid with combinations such that -

1. Digits that feature in each column or row do not appear more than once.
2. Letters that feature in each column or row do not appear more than once.
3. Each letter-digit combination does not appear more than once in the whole grid.

WARNING: Spoilers ahead, so if you want to try and solve this on your own, you might want to stop reading right now and get back after you’ve attempted it.

## The obvious and the easy

It is obvious that $$N = 2$$ has no solution. $$N = 3$$ is rather easy -

$$\begin{array}{ccc} A1 & B2 & C3 \\ B3 & C1 & A2 \\ C2 & A3 & B1 \end{array}$$

With that out of the way, now try $$N = 4$$. Really, go and try it out now. If it feels strange that the tactic for $$N = 3$$ doesn’t seem to work for $$N = 4$$, try alternatives. Yes, there are solutions for $$N = 4$$ so your attempts won’t be fruitless.

One thing to notice is that if you have one solution, you can generate other solutions by -

1. Exchanging any two columns
2. Exchanging any two rows
3. Permuting the letters across the whole board. Ex: substitute $$ABC\to BCA$$.
4. Permuting the digits across the whole board. Ex: substitute $$123\to 231$$.

Due to the above, it is easy to see that if a solution exists, you can manipulate it such that the first row will read $$A1\ B2\ C3\ldots$$ and the first coumn’s letters are also in the order $$A\ B\ C\ldots$$ when read from top to bottom. This is because you can permute the letters in such a way that the first row letters read in order, likewise for the digits and then swap rows in the appropriate order so that the left column letters read in order.

## Modulo arithmetic

At first, I tried to produce the required permutation series using modulo arithmetic … using the Julia function -

grid(N,a,b,c,d) = 1 .+ [(mod(a*i+b*j,N), mod(c*i+d*j,N)) for i in 1:N, j in 1:N]

In the above program, mod(m,n) is the remainder when you divide m by n, even considering negative values for m. For example, mod(5,3) is 2 because $$5=3+2$$, but mod(-5,3) is 1 since $$-5 = -6 + 1$$. Note that we’ll get a matrix whether we range over 1:N or 0:(N-1) for i and j.

If you try that out, you’ll find that grid(N, 1, 1, 1, 2) produces a solution for all odd values of $$N$$. However, no values of a,b,c,d work to solve for even $$N$$. You might’ve noticed that when you tried to solve for $$N = 4$$ by hand.

In the above program, I’m representing the entry for each cell as a pair of numbers (i,j) where i and j are in the range $$[1,N]$$. It is a trivial matter to translate that to a matrix of letter-digit combinations with a function like below -

function printgrid(grid)
n = size(grid)[1]
letters = "ABCDEFGHI"
digits = "123456789"
ld(i,j) = letters[grid[i,j][1]] * digits[grid[i,j][2]]
for i in 1:n
println(join([ld(i,j) for j in 1:n], " "))
end
end

## Do we have a solution?

It’s good to have a test function on a grid that we can use to check whether the grid solves the puzzle or not. To do that, we code a function that produces a boolean value of true or false depending on whether the grid solves the puzzle or not.

# We assume that the given grid is NxN in structure and
# has tuples with values in the range 0:(N-1). We won't test
# for those conditions for simplicity.
function solves(N, grid)
colsokL = all(allunique(grid[i,j][1] for j in 1:N) for i in 1:N)
colsokD = all(allunique(grid[i,j][2] for j in 1:N) for i in 1:N)
rowsokL = all(allunique(grid[j,i][1] for j in 1:N) for i in 1:N)
rowsokD = all(allunique(grid[j,i][2] for j in 1:N) for i in 1:N)
colsokL && colsokD && rowsokL && rowsokD && allunique(grid)
end

We can now use the solves function to test whether the previous grid function produces solutions for various values of N,a,b,c,d.

julia> solves(5, grid(5,1,1,1,2))
true

julia> solves(4, grid(4,1,1,1,2))
false

julia> solves(5, grid(5,1,1,1,1))
false

julia> printgrid(grid(5,1,1,1,2))
C4 D1 E3 A5 B2
D5 E2 A4 B1 C3
E1 A3 B5 C2 D4
A2 B4 C1 D3 E5
B3 C5 D2 E4 A1

Play around with the values for a,b,c,d to see what works and what doesn’t.

## Finding solutions by searching

We saw that the 4x4 case wasn’t solvable using our modulo arithmetic grid function. So how do we search for solutions in such cases?

First off, it is useful to narrow down possibilities since otherwise we’ll have to look through $$16! = 20922789888000$$ cases. We’ll use the “standard form” we came up with earlier for that - i.e. where -

1. The first row is assumed to read $$[A1,B2,C3,\ldots]$$, and
2. the first column’s letters have to be in order $$[A,B,C,\ldots]$$ when read top-down.

We’ll try and fill our the grid from top left corner and move right and down. So we’ll need some way to enumerate the possibilities at any given position assuming that all previous positions are filled in a consistent manner.

# We assume that all the elements up to (i,j) are consistent
# with the constraints.
function satisfies(N, grid, i, j, a, b)
colokL = a ∉ (grid[r,j][1] for r in 1:i-1)
colokD = b ∉ (grid[r,j][2] for r in 1:i-1)
rowokL = a ∉ (grid[i,c][1] for c in 1:j-1)
rowokD = b ∉ (grid[i,c][2] for c in 1:j-1)
ix = i*N+j
lduniq = (a,b) ∉ (grid[r,c] for r in 1:i for c in 1:N if r*N+c < i*N+j)
return colokL && colokD && rowokL && rowokD && lduniq
end

function possibilities(N, grid, i, j)
if i == 1
return [(j,j)]
end

[(r,c) for r in 1:N for c in 1:N
if (j > 1 || r == i) && satisfies(N, grid, i, j, r, c)]
end

Now we just need to explore the possibilities one by one and produce a solution if one exists.

function explore(N, grid, i, j)
if i > N return true end

for rc in possibilities(N, grid, i, j)
grid[i,j] = rc
jnext = mod(j,N)+1
inext = i + (if (jnext == 1) 1 else 0 end)
if explore(N, grid, inext, jnext)
return true
end
end
return false
end

If the above explore function returns true, then the grid matrix contains a solution. If it returns false, then there is no solution. Let’s try them in a few known situations -

julia> g = [(0,0) for i in 1:3, j in 1:3]
3×3 Matrix{Tuple{Int64, Int64}}:
(0, 0)  (0, 0)  (0, 0)
(0, 0)  (0, 0)  (0, 0)
(0, 0)  (0, 0)  (0, 0)

julia> explore(3, g, 1, 1)
true

julia> printgrid(g)
A1 B2 C3
B3 C1 A2
C2 A3 B1

julia> g = [(0,0) for i in 1:4, j in 1:4]
4×4 Matrix{Tuple{Int64, Int64}}:
(0, 0)  (0, 0)  (0, 0)  (0, 0)
(0, 0)  (0, 0)  (0, 0)  (0, 0)
(0, 0)  (0, 0)  (0, 0)  (0, 0)
(0, 0)  (0, 0)  (0, 0)  (0, 0)

julia> explore(4, g, 1, 1)
true

julia> printgrid(g)
A1 B2 C3 D4
B3 A4 D1 C2
C4 D3 A2 B1
D2 C1 B4 A3

So yeah! It was able to find the solution for the 4x4 grid.

## The biggie!

Try to solve the 6x6 grid by hand. Go ahead, really give it a try. When I did it, it felt like this problem was one of those life consuming problems like the Collatz conjecture. I later learnt (from my kid’s math teacher) that Euler struggled with the 6x6 too, so if you’re struggling, don’t be disheartened.

… or if you’re like me, you can run the explore function on it … just that you’ll have to wait a bit to get the solution.

julia> g = [(0,0) for i in 1:6, j in 1:6]
6×6 Matrix{Tuple{Int64, Int64}}:
(0, 0)  (0, 0)  (0, 0)  (0, 0)  (0, 0)  (0, 0)
(0, 0)  (0, 0)  (0, 0)  (0, 0)  (0, 0)  (0, 0)
(0, 0)  (0, 0)  (0, 0)  (0, 0)  (0, 0)  (0, 0)
(0, 0)  (0, 0)  (0, 0)  (0, 0)  (0, 0)  (0, 0)
(0, 0)  (0, 0)  (0, 0)  (0, 0)  (0, 0)  (0, 0)
(0, 0)  (0, 0)  (0, 0)  (0, 0)  (0, 0)  (0, 0)

julia> @time explore(6, g, 1, 1)
42.486919 seconds (62.89 M allocations: 5.005 GiB, 0.31% gc time)
false

What? Oh! The 6x6 has no solutions! … and we have a computational proof for that too? Apparently this took Euler ages to conjecture, which is how it started for me initially – i.e. trying to get $$6\times 6$$ to work by hand – but then I increasingly suspected there is no solution based on the patterns I was seeing, then decided to write a program … etc.

Note that we weren’t particularly aiming for efficient code, just obvious code, so I expect we can make this faster as well. Maybe that’s the topic of another post, but you can get a quick win by just resorting to short-circuit evaluation for the satisfies function -

function satisfies(N, grid, i, j, a, b)
ix = i*N+j
(a ∉ (grid[r,j][1] for r in 1:i-1)
&& b ∉ (grid[r,j][2] for r in 1:i-1)
&& a ∉ (grid[i,c][1] for c in 1:j-1)
&& b ∉ (grid[i,c][2] for c in 1:j-1)
&& (a,b) ∉ (grid[r,c] for r in 1:i for c in 1:N if r*N+c < i*N+j))
end

The idea is that if the first condition is false, then the subsequent ones don’t need to be checked, thereby saving time. This version, on my machine runs in < 10 seconds - i.e. gives about a 4x speed up already. A larger thought worth reflecting on is that what would’ve taken Euler perhaps a lifetime to check can now be checked in under 10 seconds! We have a computer proof on our hands … as long as we’re confident enough that our program is correct.

## The 4x4 pattern

So now we have an easy solution for odd $$N$$, and a weird looking solution for $$4\times 4$$ and we also know that there is no solution for $$6\times 6$$. Is there a solution for $$8\times 8$$?

When considering this, my thinking took a drastic turn. I tried to extract a pattern out of the 4x4 solution - a single instance - and generalize it for the 8x8. What was odd was that it worked!

Let’s look at the 4x4 solution again -

$$\begin{array}{cccc} A1 & B2 & C3 & D4 \\ B3 & A4 & D1 & C2 \\ C4 & D3 & A2 & B1 \\ D2 & C1 & B4 & A3 \end{array}$$

If we focus on the letter patterns, we see the permutations - $$ABCD$$, $$BADC$$, $$CDBA$$ and $$DCBA$$. These looked very familiar to me. They reminded me of … the Walsh-Hadamard transform matrix –

$$\begin{array}{rrrr} 1 & 1 & 1 & 1 \\ 1 & -1 & 1 & -1 \\ 1 & 1 & -1 & -1 \\ 1 & -1 & -1 & 1 \end{array}$$

… but that matrix didn’t seem to have anything to do with permutations. Ok let’s look at the digit patterns, which read – $$1234$$, $$3412$$, $$4321$$ and $$2143$$. This is essentially the same permutation sequence as with the letters, but in a different order. If we number the permutation sequence of the letters as $$0123$$, this order is $$0231$$.

So, what does the Walsh-Hadamard matrix for an 8x8 transform look like?

$$\begin{array}{rrrrrrrr} 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 & -1 & -1 & -1 & -1 & 1 & 1 \\ 1 & -1 & -1 & 1 & -1 & 1 & 1 & -1 \end{array}$$

How to make sense of this as a set of permutations? One thing that can help is to note that the rows of the matrix can be generated from the following fundamental patterns -

$$\begin{array}{lcrrrrrrrr} P1 & = & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ P2 & = & 1 & 1 & 1 & 1 & -1 & -1 & -1 & -1 \\ P3 & = & 1 & 1 & -1 & -1 & 1 & 1 & -1 & -1 \\ P4 & = & 1 & -1 & 1 & -1 & 1 & -1 & 1 & -1 \end{array}$$

Given that we’d take the first row there - $$P1$$ - to mean “no permutations”, we can set that aside for the moment.

The rows of the 8x8 matrix can be expressed in terms of these as -

$$\begin{array}{l} P1 \\ P4 \\ P3 \\ P3 \times P4\ \text{(element-wise multiplication)} \\ P2 \\ P2 \times P4 \\ P2 \times P3 \\ P2 \times P3 \times P4 \end{array}$$

Note that the second half is the same as the first half, except with $$P2 \times$$ as an added permutation. If we take the four permutations $$P1,P2,P3,P4$$ to mean $$ABCDEFGH$$, $$EFGHABCD$$, $$CDABGHEF$$ and $$BADCFEHG$$ respectively, we can use the multiplication to mean application of permutations (which in this case are also commutative) to get the following letter sequence –

$$\begin{array}{lcccccccc} 0) & A & B & C & D & E & F & G & H \\ 1) & B & A & D & C & F & E & H & G \\ 2) & C & D & A & B & G & H & E & F \\ 3) & D & C & B & A & H & G & F & E \\ 4) & E & F & G & H & A & B & C & D \\ 5) & F & E & H & G & B & A & D & C \\ 6) & G & H & E & F & C & D & A & B \\ 7) & H & G & F & E & D & C & B & A \end{array}$$

The above pattern meets the “no repeated letters in columns and rows” criterion. We need to pick a different ordering of this pattern for the digits. With some fooling around, I tried $$04627315$$, which gave me –

$$\begin{array}{cccccccc} A1 & B2 & C3 & D4 & E5 & F6 & G7 & H8 \\ B5 & A6 & D7 & C8 & F1 & E2 & H3 & G4 \\ C7 & D8 & A5 & B6 & G3 & H4 & E1 & F2 \\ D3 & C4 & B1 & A2 & H7 & G8 & F5 & E6 \\ E8 & F7 & G6 & H5 & A4 & B3 & C2 & D1 \\ F4 & E3 & H2 & G1 & B8 & A7 & D6 & C5 \\ G2 & H1 & E4 & F3 & C6 & D5 & A8 & B7 \\ H6 & G5 & F8 & E7 & D2 & C1 & B4 & A3 \end{array}$$

.. which solves the 8x8 case. Now what exactly is the $$04627315$$ pattern? If we write those numbers out in binary (3 bits), we get –

$$\begin{array}{c} 000 \\ 100 \\ 110 \\ 010 \\ \\ 111 \\ 011 \\ 001 \\ 101 \end{array}$$

We see that the first four and the last four in the above sequence are bit-wise inverses. Within each set of four as well, if you drop the last bit, you can again divide it into two where the first half and the second half are bit-wise flipped – i.e. $$00, 10$$ and $$11, 01$$ are bitwise flipped. So we have an algorithm to construct the sequence for the digits given the sequence for the letters for $$N = 2^k$$ (for which the Walsh-Hadamard matrix is defined). Let’s write it out for $$16\times 16$$, for example –

$$\begin{array}{ccc} 0000 & = & 0 \\ 1000 & = & 8 \\ 1100 & = & 12 \\ 0100 & = & 4 \\ 1110 & = & 14 \\ 0110 & = & 6 \\ 0010 & = & 2 \\ 1010 & = & 10 \\ 1111 & = & 15 \\ 0111 & = & 7 \\ 0011 & = & 3 \\ 1011 & = & 11 \\ 0001 & = & 1 \\ 1001 & = & 9 \\ 1101 & = & 13 \\ 0101 & = & 5 \end{array}$$

Just to be sure, we need to check whether this indeed solves the $$16\times 16$$ case. So let’s work that out –

A1  B2  C3  D4  E5  F6  G7  H8  I9  J10 K11 L12 M13 N14 O15 P16
B9  A10 D11 C12 F13 E14 H15 G16 J1  I2  L3  K4  N5  M6  P7  O8
C13 D14 A15 B16 G7  H8  E9  F10 K5  L6  I7  J8  O1  P2  M3  N4
D5  C6  B7  A8  H1  G2  F3  E4  L13 K14 J15 I16 P7  O8  N9  M10
E15 F16 G13 H14 A9  B10 C11 D12 M7  N8  O5  P6  K3  L4  I1  J2
F7  E8  H5  G6  B3  A4  D1  C2  N15 M16 P13 O14 L9  K10 J11 I12
G3  H4  E1  F2  C7  D8  A5  B6  O11 P12 M9  N10 I15 J16 K13 L14
H11 G12 F9  E10 D15 C16 B13 A14 P3  O4  N1  M2  J7  I8  L5  K6
I16 J15 K14 L13 M10 N9  O12 P11 A8  B7  C6  D5  E4  F3  G2  H1
J8  I7  L6  K5  N4  M3  P2  O1  B16 A15 D14 C13 F12 E11 H10 G9
K4  L3  I2  J1  O8  P7  M6  N5  C12 D11 A10 B9  G16 H15 E14 F13
L12 K11 J9  I10 P16 O15 N14 M13 D4  C3  B2  A1  H8  G7  F6  E5
M2  N1  O4  P3  K8  L7  I6  J5  E10 F9  G12 H11 A14 B13 C16 D15
N10 M9  P12 O11 L14 K13 J16 I15 F2  E1  H4  G3  B6  A5  D8  C7
O14 P13 M16 N15 I12 J11 K10 L9  G6  H5  E8  F7  C2  D1  A4  B3
P6  O5  N8  M7  J2  I1  L4  K3  H13 G13 F16 E15 D12 C11 B10 A9

Well, that looks fine to me, but I admit I may miss repetitions there and so we need to really check it with code … which I’ll set aside for now.1 2

So now, how do we compute that permutation sequence for the digits – i.e. produce the index table above.

function dmap(N,i)
if N == 2 return i end

N2 = div(N,2)
if i < N2
dmap(N2, i) * 2
else
xor(N2-1, dmap(N2,i-N2)) * 2 + 1
end
end

Now we can test whether we get the same sequence for 16 –

julia> dmap.(16,0:15)
16-element Vector{Int64}:
0
8
12
4
14
6
2
10
15
7
3
11
1
9
13
5

So there you go. It looks like we have a way to get a solution for the case where $$N = 2^k$$.

## Proofs, Prooofs, Proooofs!

What we have here is a bunch of conjectures. How do we know that this scheme for $$N = 2^k$$ will work?

If we consider how we built up the permutations from the Hadamard matrix, we can rewrite it rather systematically as below (for $$8\times 8$$) –

$$\begin{array}{lccc} \text{permutation} & \text{bit pattern} & \text{decimal} & \text{dmap}(8,i) \\ P2^0 \times{} P3^0 \times{} P4^0 & 000 & 0 & 0 \\ P2^0 \times{} P3^0 \times{} P4^1 & 001 & 1 & 4 \\ P2^0 \times{} P3^1 \times{} P4^0 & 010 & 2 & 6 \\ P2^0 \times{} P3^1 \times{} P4^1 & 011 & 3 & 2 \\ P2^1 \times{} P3^0 \times{} P4^0 & 100 & 4 & 7 \\ P2^1 \times{} P3^0 \times{} P4^1 & 101 & 5 & 3 \\ P2^1 \times{} P3^1 \times{} P4^0 & 110 & 6 & 1 \\ P2^1 \times{} P3^1 \times{} P4^1 & 111 & 7 & 5 \end{array}$$

So the permutation set for the letters is straight forward, and since we know that the permutations $$P2$$, $$P3$$ and $$P4$$ commute and are also orthogonal, we can see that the above set will produce no overlaps for the letters when used to make a matrix.

So, given a sequence of length $$N = 2^k$$, how do we produce the “Walsh-Hadamard permutation sequence” for it? To start with, we need to compute the “basis permutations” –

# k ranges over 2^k < N and N = 2^m
function pk(N, k)
p = 2^k
p2 = 2*p
if p2 == N
return vcat(p+1:N, 1:p)
end
[i - mod(i,p2) + pk(p2,k)[mod(i,p2)+1] for i in 0:N-1]
end

# N must be a power of 2
perms = []
seq = vcat(1:N)
for i in 0:(N-1)
s = copy(seq)
(j,n) = (i,0)
while 2^n < N
s = if mod(j,2) == 1 s[pk(N,n)] else s end
(j,n) = (div(j,2), n+1)
end
push!(perms, s)
end
perms
end

Julia hint: Each entry of the result of walsh_hadamard_permutations will be an index array that you can use to index a sequence to permute it.

That was a bit circuitous actually, as the permutations can be generated more easily recursively like this, where we’re producing a matrix instead of an array of arrays -

function walsh_hadamard_permutations(N)
if N == 1
return [1]
end

N2 = div(N,2)
vcat(hcat(p, p .+ N2), hcat(p .+ N2, p))
end

Now, why would using the dmap approach produce a digit pattern such that we get no repeats for the letter-digit combinations? I’m not sure I know how to prove that for the general $$N = 2^k$$ case yet. Maybe this is indeed a life spoiler problem and not knowing the reason will continue to haunt me for a while. For now, I’m settling for “it looks reasonable and within reach” and am satisfied with the fact that I found a reordering sequence for $$[0,2^k)$$ that I didn’t know about before this, purely by playing around with patterns.

## (21 Feb 2022) Continuing the journey

In 1 and 2, I noted after writing the material above that the 16x16 I had posted there was wrong and wasn’t a solution and also conjectured that there are no solutions for $$N = 2^{3k+1}$$. Well, I explored the Walsh-Hadamard based technique some more and after one more generalization step was able to produce (and verify) a solution for the 16x16 which I give below without much fanfare (so much for that conjecture!) –

A1  B2  C3  D4  E5  F6  G7  H8  I9  J10 K11 L12 M13 N14 O15 P16
B5  C6  D7  A8  F9  G10 H11 E12 J13 K14 L15 I16 N1  O2  P3  M4
C9  D10 A11 B12 G13 H14 E15 F16 K1  L2  I3  J4  O5  P6  M7  N8
D13 A14 B15 C16 H1  E2  F3  G4  L5  I6  J7  K8  P9  M10 N11 O12
E6  F7  G8  H5  I10 J11 K12 L9  M14 N15 O16 P13 A2  B3  C4  D1
F10 G11 H12 E9  J14 K15 L16 I13 N2  O3  P4  M1  B6  C7  D8  A5
G14 H15 E16 F13 K2  L3  I4  J1  O6  P7  M8  N5  C10 D11 A12 B9
H2  E3  F4  G1  L6  I7  J8  K5  P10 M11 N12 O9  D14 A15 B16 C13
I11 J12 K9  L10 M15 N16 O13 P14 A3  B4  C1  D2  E7  F8  G5  H6
J15 K16 L13 I14 N3  O4  P1  M2  B7  C8  D5  A6  F11 G12 H9  E10
K3  L4  I1  J2  O7  P8  M5  N6  C11 D12 A9  B10 G15 H16 E13 F14
L7  I8  J5  K6  P11 M12 N9  O10 D15 A16 B13 C14 H3  E4  F1  G2
M16 N13 O14 P15 A4  B1  C2  D3  E8  F5  G6  H7  I12 J9  K10 L11
N4  O1  P2  M3  B8  C5  D6  A7  F12 G9  H10 E11 J16 K13 L14 I15
O8  P5  M6  N7  C12 D9  A10 B11 G16 H13 E14 F15 K4  L1  I2  J3
P12 M9  N10 O11 D16 A13 B14 C15 H4  E1  F2  G3  L8  I5  J6  K7

I’ll take a moment there to share the combined thrill of finding a solution and the disappointment of having made a failed conjecture. Well, apparently it’s been proved that solutions exist for all N except for $$N = 2\ \text{and} 6$$!! Isn’t that weird? For the moment I’m going to ignore that and pretend that there is no such proof available .. at least not one that proves by providing a construction.

So how did that solution come about??

### A generalization of the “Walsh-Hadamard permutations”

I wanted to see whether the “Walsh-Hadamard permutations” can be extended to non-powers of 2. To get there, let’s look at the power-of-2 one I’d given earlier -

function walsh_hadamard_permutations(N)
if N == 1
return [1]
end

N2 = div(N,2)
vcat(hcat(p, p .+ N2), hcat(p .+ N2, p))
end

If we look at what’s happening in the recursive step, we see a division by 2 and a replication of the “N/2” permutations in the last line in a 2x2 pattern. What if we parameterize this “2” and control what is used in every recursive descent step? To get at that, we can change the definition of the parameter $$N$$ to be a “vector of bases”, each entry indicating what value to use at each successive recursive descent step. Also, it would be nice to be able to calculate single entries of the resultant permutation matrix, or single rows, without having to calculate the entire matrix. So I remodelled the function to produce a function of two indices (a “row index” and a “column index”) which when called with the appropriate indices yields the matrix’s value.

For non-2 values, I used a generalization of the bit swap operation “1,2;2,1” that you get with a conventional Walsh-Hadamard, to be a successive rotation - like “1,2,3;2,3,1;3,1,2” for example, when we produce a 3x3 permutation matrix. So putting all of that together, I generalized the above approach to yield -

# "gwhp = generalized Walsh-Hadamard permutations"
function gwhp(N)
if length(N) == 0
return (i,j) -> 0
# Returning 1 in this function will produce a matrix whose entries
# for a 1-based index instead of a 0-based index.
end

bn = N[end]
N2 = N[1:end-1]
bn2 = prod(N2)
wh2 = gwhp(N2) # Calcuating it here instead of inside the function
# will avoid the recursive step when calculating the
# value of each matrix entry.

(r,c) -> wh2(mod(r,bn2), mod(c,bn2)) + bn2 * mod(div(c,bn2) + div(r,bn2), bn)
# The first part of the expression produces the same pattern for
# each block of bn2 x bn2. The second part varies each block according
# to the rotation pattern. Note that the values of r and c in the second
# part will take on constant values for each bn2 x bn2 block. Together,
# they achieve the same effect as the vcat and hcat in the earlier
end

# We need to take n as a separate parameter since when given
# function representation we can't tell what its intended size is.
function printwhp(whp, n)
for i in 0:n-1
for j in 0:n-1
print(whp(i,j), " ")
end
println()
end
end

One additional thing to note is that the indices are zero-based and not one-based … since we’re using modulo arithmetic here and zero-based is easier to work with when using modulo arithmetic. The entries of the permutation matrix are also therefore zero-based. Let’s see how to get the 8x8 Walsh-Hadamard permutations matrix -

julia> printwhp(gwhp([2,2,2]), 8)
0 1 2 3 4 5 6 7
1 0 3 2 5 4 7 6
2 3 0 1 6 7 4 5
3 2 1 0 7 6 5 4
4 5 6 7 0 1 2 3
5 4 7 6 1 0 3 2
6 7 4 5 2 3 0 1
7 6 5 4 3 2 1 0

To see what it does for odd values, here are a few examples –

julia> printwhp(gwhp([5]), 5)
0 1 2 3 4
1 2 3 4 0
2 3 4 0 1
3 4 0 1 2
4 0 1 2 3

julia> printwhp(gwhp([2,3]), 6)
0 1   2 3   4 5
1 0   3 2   5 4

2 3   4 5   0 1
3 2   5 4   1 0

4 5   0 1   2 3
5 4   1 0   3 2

julia> printwhp(gwhp([3,2]), 6)
0 1 2   3 4 5
1 2 0   4 5 3
2 0 1   5 3 4

3 4 5   0 1 2
4 5 3   1 2 0
5 3 4   2 0 1

Note: I’ve added extra spaces to make the patterns stand out in the above examples. printwhp itself won’t print those spaces. You can try to modify it to do that if you’re interested.

So having generalized the WHPs thus, we also need to generalize some of the other functions before we can get a letter-digit matrix computed for non-power-of-2 values of N – dmap, and a function that puts together gdmap (the generalization of dmap) and gwhp to produce the result matrix.

The generalization of dmap is interesting. We need to find a generalization for the “bit flip” operation that is recursively applied. To do this, we treat the vector $$N$$ as the “bases” of a “generalized place value system” $$N = [b_1,b_2,b_3,\ldots]$$ where we express any number $$n$$ as a unique series of “generalized digits” $$[d_1,d_2,d_3,\ldots]$$ such that $$0 \leq d_k \lt b_k$$ and $$n = \sum_{k}{d_k \prod_{0}^{k-1}{b_m}}$$ where we take $$b_0 = 1$$. If all our $$b_k = 2$$, we can see that the “bit flip” operation is a digit-wise increment $$\mod b_k$$. So that’s what we’ll take as the generalization of the “bit flip”. This gives us the following “generalized dmap” –

"""
brot(N, m, n, f=1)

N is a "generalized number base" vector.
m is the number of rotation steps to applu to each digit
(modulo the corresponding base).
n is the number whose digits are to be "rotated".

If N is a power of 2 and $fm != 0 \mod 2$, then this is
effectively just xor(N-1, n) which flips all the bits
in the binary representation of n.
"""
function brot(N, m, n, f=1)
if m == 0 return n end
rn = 0
bs = 1
for bk in N
dk = mod(n, bk)

# rn rebuilds the number n but with the digits
# "rotated" by m steps
rn += bs * mod(dk+f*m,bk)

n = div(n, bk)
bs *= bk
end
rn
end

function gdmap(N, f=1)
if length(N) == 1 return (i -> mod(f*i,N[1])) end

N2 = N[1:end-1]
bn2 = prod(N2)
g2 = gdmap(N2,f)

function mi(i)
m = div(i, bn2)
brot(N2, m, g2(i - m * bn2)) * N[end] + m
end
end

We’ve also used the same approach of producing a function instead of a vector, so we can compute individual offsets separately without having to compute potentially large vectors and matrices.

Correspondingly, we can define the generalized version of the function that combines the gwhp and gdmap to produce the letter-digit matrix as –

function gmat(N,f=1)
w = gwhp(N)
d = gdmap(N,f)

# The digits are the same pattern as the letters except for a reordering
# of the rows.
(r,c) -> (w(r,c), w(d(r),c))
end

Some examples of what gmat produces –

julia> printgrid(gmat([2,2]),4)
A1 B2 C3 D4
B3 A4 D1 C2
C4 D3 A2 B1
D2 C1 B4 A3

julia> printgrid(gmat([2,3]),6)
A1 B2 C3 D4 E5 F6
B4 A3 D6 C5 F2 E1
C5 D6 E1 F2 A3 B4
D2 C1 F4 E3 B6 A5
E3 F4 A5 B6 C1 D2
F6 E5 B2 A1 D4 C3

julia> printgrid(gmat([3,2]),6)
A1 B2 C3 D4 E5 F6
B3 C1 A2 E6 F4 D5
C5 A6 B4 F2 D3 E1
D6 E4 F5 A3 B1 C2
E2 F3 D1 B5 C6 A4
F4 D5 E6 C1 A2 B3

… where we adapt the previous printgrid for the new function representation like this -

function printgrid(m, n)
letters = "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
digits = string.(1:26)
ld(i,j) = letters[1+m(i,j)[1]] * digits[1+m(i,j)[2]]
for i in 0:n-1
println(join([ld(i,j) for j in 0:n-1], " "))
end
end

So the gmat function produces the following for gmat([4,4])

julia> printgrid(gmat([4,4], -1),16)
# .... what I showed at the start of this section ...

julia> printgrid(gmat([4,4]), 16)
A1  B2  C3  D4  E5  F6  G7  H8  I9  J10 K11 L12 M13 N14 O15 P16
B5  C6  D7  A8  F9  G10 H11 E12 J13 K14 L15 I16 N1  O2  P3  M4
C9  D10 A11 B12 G13 H14 E15 F16 K1  L2  I3  J4  O5  P6  M7  N8
D13 A14 B15 C16 H1  E2  F3  G4  L5  I6  J7  K8  P9  M10 N11 O12
E14 F15 G16 H13 I2  J3  K4  L1  M6  N7  O8  P5  A10 B11 C12 D9
F2  G3  H4  E1  J6  K7  L8  I5  N10 O11 P12 M9  B14 C15 D16 A13
G6  H7  E8  F5  K10 L11 I12 J9  O14 P15 M16 N13 C2  D3  A4  B1
H10 E11 F12 G9  L14 I15 J16 K13 P2  M3  N4  O1  D6  A7  B8  C5
I11 J12 K9 L10  M15 N16 O13 P14 A3  B4  C1  D2  E7  F8  G5  H6
J15 K16 L13 I14 N3  O4  P1  M2  B7  C8  D5  A6  F11 G12 H9  E10
K3  L4  I1  J2  O7  P8  M5  N6  C11 D12 A9  B10 G15 H16 E13 F14
L7  I8  J5  K6  P11 M12 N9  O10 D15 A16 B13 C14 H3  E4  F1  G2
M8  N5  O6  P7  A12 B9  C10 D11 E16 F13 G14 H15 I4  J1  K2  L3
N12 O9  P10 M11 B16 C13 D14 A15 F4  G1  H2  E3  J8  K5  L6  I7
O16 P13 M14 N15 C4  D1  A2  B3  G8  H5  E6  F7  K12 L9  I10 J11
P4  M1  N2  O3  D8  A5  B6  C7  H12 E9  F10 G11 L16 I13 J14 K15

The second 16x16 also satisfies the constraints and is a solution.

We can also adapt the solves function for the new representation to use less memory and rely more on compute instead, as follows –

function solves(m, N)
ixs = 0:N-1
row(r) = [m(r,c) for c in ixs]
(all(allunique(m(i,j)[1] for j in ixs) for i in ixs)
&& all(allunique(m(i,j)[2] for j in ixs) for i in ixs)
&& all(allunique(m(j,i)[1] for j in ixs) for i in ixs)
&& all(allunique(m(j,i)[2] for j in ixs) for i in ixs)
&& all(allunique(vcat(row(r1),row(r2))) for r1 in ixs for r2 in ixs if r2 < r1))
end

Using the above solves, we can check whether we have a solution –

julia> solves(gmat([4,4]), 16)
true

julia> solves(gmat([2,3]), 6)
false

julia> solves(gmat([2,2,2,2]), 16)
false

julia> solves(gmat([2,2,2]), 8)
true

julia> solves(gmat([3,5]), 15)
true

julia> solves(gmat([5]), 5)
false

The last result is rather disappointing. the matrix produced reads -

A1 B2 C3 D4 E5
B2 C3 D4 E5 A1
C3 D4 E5 A1 B2
D4 E5 A1 B2 C3
E5 A1 B2 C3 D4

Aha - while our brot corresponds to “bit flip” when the bases are 2, it doesn’t quite flip anything for other bases. All we have to do to get a solution is to pass in an f value of -1 … which produces –

julia> printgrid(gmat([5],-1),5)
A1 B2 C3 D4 E5
B5 C1 D2 E3 A4
C4 D5 E1 A2 B3
D3 E4 A5 B1 C2
E2 A3 B4 C5 D1

julia> solves(gmat([5],-1),5)
true

So we now have a generalized approach that solves more than just the odd sized matrix.

I’m getting more interested in proofs after this much exploration and explanation … so the journey into patterns shall continue I suppose!

## (22-02-2022) Proving the “odd N” case

Let’s start with the simple case - of odd $$N$$, i.e. $$N \equiv 1 \mod 2$$. We constructed the matrix using the function -

$$\text{grid}(i,j) = ((i+j\mod N), (i+2j\mod N))$$

The first and second components generate the whole range $$[0,N)$$ only for odd $$N$$. So for even $$N$$, the rows and/or columns will have repeated values due to the stepping by 2. Or to put it in modulo arithmetic terms,

$$\begin{array}{lcccr} & 2i & \equiv & 2j & \mod N \\ \implies & i & \equiv & j & \mod N \\ \text{iff} & \text{gcd}(2,N) & = & 1 \end{array}$$

So what we need to find out is whether $$\text{grid}(i,j) = \text{grid}(m,n)$$ implies $$i = m$$ and $$j = n$$.

The equation gives us –

$$\begin{array}{llcr} & i+j & \equiv & m+n\mod N \\ \text{and} & i+2j & \equiv & m+2n\mod N \\ \implies\text{ (by subtraction)} & j & \equiv & n\mod N \end{array}$$

And since $$0 \leq j \lt N$$ and $$0 \leq n \lt N$$, we have $$j = n$$. Similarly, we can show that $$i = m$$. So we now know why $$\text{grid}(i,j)$$ produces a solution.

You might’ve noticed that we skipped the row and column uniqueness criteria for letters and digits. Those are trivially ensured by the function since for a given $$i$$, the columns all take on different values and likewise for a given $$j$$, the rows all take on different values – for both components.

### Why can’t we apply that for the “even N” case?

From the previous argument for “odd N”, it looks like we need to find appropriate $$a,b,c,d$$ that meet the criteria for even $$N$$. For example, if we take $$N = 8$$, we have two numbers $$3$$ and $$5$$ which are both coprime to $$8$$. So what do we get if we construct a $$\text{grid}(i,j)$$ function based on these two numbers like shown below –

$$\text{grid}_8(i,j) = (i + 3j\ \mod 8, i + 5j\ \mod 8)$$

With the above function, it is easy to see that we’ll be able to meet all of the row and column uniqueness criteria, since $$\text{gcd}(3,8) = 1$$ and $$\text{gcd}(5,8) = 1$$.

However, when we set out to show that $$\text{grid}_8(i,j) = \text{grid}_8(m,n) \implies i = m\ \text{and}\ j = n$$, we will need to show that from the following equations -

$$\begin{array}{rclr} i + 3j & \equiv & m + 3n & \mod 8\\ i + 5j & \equiv & m + 5n & \mod 8 \end{array}$$

When we subtract the two equations, we get $$2j \equiv 2n\ \mod 8$$. That equation does not imply $$j = n$$ because $$\text{gcd}(2,8) \neq 1$$. This means there exist two different values for $$j$$ and $$n$$ which can be used to satisfy $$2j \equiv 2n\ \mod 8$$, even with the constraints $$0 \leq j \lt 8$$ and $$0 \leq n \lt 8$$. For example, $$2\times 3 \equiv 2\times 7\ \mod 8$$ because $$2\times 7 = 14$$ which leaves a remainder $$6$$ when divided by $$8$$. Similarly, $$2\times 2 \equiv 2\times 6\ \mod 8$$ since $$12$$ leaves a remainder of $$4$$ when divided by $$8$$.

Given that $$\text{gcd}(k,N) = 1$$ where $$N$$ is even, implies that $$k$$ is odd (otherwise the GCD will be at least 2), it implies that using two odd coefficients can never yield a unique-element matrix because the difference between two odd numbers is always even which will have a GCD of at least 2 with an even number.

So what about using values other than 1 for $$a$$ and $$c$$?

Our equations then become -

$$\begin{array}{rclr} ai + bj & \equiv & am + bn & \mod N \\ ci + dj & \equiv & cm + dn & \mod N \end{array}$$

We know that $$i \equiv j\ \mod N$$ implies $$ai \equiv aj\ \mod N$$ for all $$a$$. So we can multiply the equations by opposite coefficients to get –

$$\begin{array}{rclr} cai + cbj & \equiv & cam + cbn & \mod N \\ aci + adj & \equiv & acm + adn & \mod N \end{array}$$

If we then subtract the equations, we get –

$$(ad-bc)j \equiv (ad-bc)n\ \mod N$$

For that to imply $$j \equiv n\ \mod N$$, we require $$\text{gcd}(ad-bc,N) = 1$$. That implies that both $$ad$$ and $$bc$$ cannot be odd, since the difference will not be coprime to $$N$$. One of them will have to be even. If one of them, say $$ad$$ is even, then that implies either $$a$$ is even or $$d$$ is even. In either case we will fail the column/row uniqueness critertion for even $$N$$.

Therefore it is not possible to use the modulo arithmetic approach to construct a solution for even $$N$$.

So the initial exploration we did amounts to treating computation as “experimental mathematics” :)

1. PS: There are errors there that my eye missed (ex: H8 in the first and 3rd rows). When I checked using the solves function, it looks like the 16 case isn’t solved by this approach, though the 8 case is. But all is not lost. After all, we know that it doesn’t solve the case of $$2^1$$ either, for which there is no solution. When I checked for a few powers of 2 up to $$2^{12}$$, I found that the approach doesn’t work for all powers of the form $$2^{4+3k}$$ for integer $$k$$ and (to my relief) works for other powers of 2, at least up to $$2^{12}$$. This gets more and more interesting! ↩

2. PS2: Given that $$2^1$$ is also of the same pattern and the approach certainly does not work for it and there is no solution for it, I’m at this point conjecturing that there are no solutions for $$N = 2^{3k+1}$$. Why those? ↩


text:
# Tag Info

2

Fixed point processing is very difficult. Floating point is A LOT easier. The best algorithm and scaling approach depends a lot on your specific filter and the statistics of your signal. There is no "one size fits all" solution. Cascaded second order sections are almost always the way to go. Primarily they guarantee that your coefficients values are bounded ...

1

I did not read carefully through the whole question as I don't have time now, but have you tried some form of robust peak detection? See e.g., https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.find_peaks.html Then you can set parameters such as the minimum (and maximum) distance, the prominence, the minimal height. If you look for minima ...

0

Hinted by the answer from @chipaudette I made a Python implementation using the s-domain transfer function (https://en.wikipedia.org/wiki/A-weighting). # python 3.7 import scipy.signal as signal import math as math import matplotlib.pyplot as plt import array as arr # for discretization sr = 48000.0 # enter the zeros, poles and scaling as in s-domain ...

3

First of all, it's not correct to say "poles should (always) be inside the unit circle for an LTI system to be stable" ; unless it's implied that system is also causal. Otherwise, if the system is noncausal, then its poles should be outside of unit circle for the system being stable. For IIR systems that are described by LCCDEs causality must be externally ...

0

From the difference equation alone you cannot tell if a system is causal or not. For example, the difference equation $$y[n]=-a_1y[n-1]-a_2y[n-2]+b_0x[n]\tag{1}$$ can be used to describe three different systems. The first one is a causal system, as suggested by the form given in $(1)$. However, rewriting $(1)$ as $$y[n-2]=-\frac{a_1}{a_2}y[n-1]-\frac{1}{... 0 Substitute m = n+K, i.e., n=m-K. This gives$$y[m-K] = -a_1 y[m-K+1] - \ldots - a_{K-1} y[m-1] - a_K y[m] + b_0 x[m-K] + b_1 x[m-K-1] + \ldots+ b_L x[m-K-L].$$Rearranging gives$$ y[m] = -\frac{a_{K-1}}{a_K} y[m-1] - \frac{a_{K-2}}{a_K} y[m-2]-\ldots-\frac{a_1}{a_K}y[m-K+1] - \frac{1}{a_K}y[m-K] + \frac{b_0}{a_K}x[m-K] + \ldots + \frac{b_L}{a_K}x[m-K-L]....

2

Much depends on what you want to deduce from your data. In general, if your sample rate is not uniform, your measurements should be accurately time stamped. The nearer to an average of 1/10 a second your sample intervals are the better. A typical hueristic used in Engineering is the 1/10 rule, so if your samples are within 1/100 of a second of 1/10 ...

1

There is no reasonable way to explicitly choose the number of taps before and after the peak. The reason is simple: the arbitrary magnitude response design results in a linear phase, and, consequently, the impulse response is symmetrical. If you specify a desired complex frequency response in terms of magnitude and phase, the location of the peak is an ...

3

In general, there is no such requirement for notch filters that $H(e^{j0})=H(e^{j\pi})$ must be satisfied. You could definitely have a notch filter with $H(e^{j0})\neq H(e^{j\pi})$. Having the same gain at DC and at Nyquist is just a practical definition, and if you have a sufficient number of degrees of freedom (i.e., filter coefficients) you might as well ...

2

Apart from the method operating on the cepstrum suggested in hotpaw2's answer, there is the obvious way of reflecting the zeros outside the unit circle into the unit circle by keeping their angle unchanged and inverting their magnitude. So each zero $z_0$ with $|z_0|>1$ needs to be replaced by $1/z_0^*$, where $*$ denotes complex conjugation. You also ...

3

I've seen Julius' MATLAB code and I know what it does. Essentially, given an LTI filter with impulse response, $h[n]$, and frequency response: \begin{align} H(e^{j \omega}) &\triangleq \Big| H(e^{j \omega}) \Big| \, e^{j \phi(\omega) } \\ &= \sum\limits_{n=-\infty}^{\infty} h[n] \, e^{-j \omega n} \end{align} Then \$\Big| H(e^{j \omega}) \Big| ...

1

mps, a matlab/octave script using Cepstral methods to estimate a minimum phase form for a spectrum is here: https://ccrma.stanford.edu/~jos/filters/Matlab_listing_mps_m_test.html#sec:tmps explanation on same (JOS) site: https://ccrma.stanford.edu/~jos/filters/Conversion_Minimum_Phase.html

Top 50 recent answers are included