text:
Codeforces celebrates 10 years! We are pleased to announce the crowdfunding-campaign. Congratulate us by the link https://codeforces.com/10years. ×

### Chilli's blog

By Chilli, history, 11 months ago, ,

In various places(12), people have talked about "Dinic's With Scaling" having $O(VElog(U))$ complexity, where U is the max edge capacity. However, in just as many places, people have cast doubt about this complexity or worried that this algorithm ends up being slower in practice(12).

I had the same doubts. In fact, a couple months ago I benchmarked Dinic's with scaling on some hard test cases and thought that Dinic's with/without scaling performed similarly.

As it turns out, I benchmarked it again yesterday and found that Dinic's with scaling definitely has significantly better asymptotic complexity (by a factor of V).

## Synthetic Benchmarking

I used the hard test case generator at 2 posted by ko_osaga here. The one posted at 1 creates a graph with only V edges and moreover, doesn't work for me... Notably, given a V, it generates a graph with V nodes and $O(V^2)$ edges.

The code I used was my own Dinic's implementation (originally from emaxx) found here. Notably, the SCALING variable toggles between Dinic's with scaling and without. It's a small change, as you can see.

I ran the two code across a variety of input sizes, averaging across 10 runs, then taking the log of input size and times, plotting them in a where blue is the runtime of Dinic's without scaling and orange is the runtime of Dinic's with scaling. Since the important thing to look at is the slope of the line, I normalized the lines to start at the same point.

Fitting a linear regression to the points, we see that the slope of Dinic's without scaling is 3.204, while the slope of Dinic's with scaling is 1.93, which correspond to complexities of $O(V^3)$ and $O(V^2)$ respectively.

In terms of actual runtime, here are some numbers on specific sizes.

## Benchmarking on Problems

There are 3 flow problems, I benchmarked on — SPOJ FASTFLOW, (VNSPOJ FASTFLOW)[https://vn.spoj.com/problems/FFLOW/), and a Chinese LOJ Flow Problem (it's intended to solve this problem with Highest Label Preflow Push, which has $O(V^2\sqrt{E})$ complexity).

Also, to vouch for my implementation, here's some benchmarks with other implementations :^) Mine (to enable scaling change the SCALING boolean), Emaxxs dinic,LOJ fastest (HLPP, the fastest submission on the Chinese flow problem). As I was doing these benchmarks I realized the LOJ submission is actually faster than mine on all of my benchmarks :'(. It is substantially longer than my implementation though :^)

SPOJ (V<=5000, E<=3e4)
Mine: 0.04
Emaxx Dinic: 0.11
Mine w/ scaling: 0.27
LOJ fastest: 0.00

VNSPOJ (V<=5000, E<=3e4, harder test cases):
Mine: TLE
Emaxx dinic: TLE
Mine w/ scaling: 0.00
LOJ fastest: 0.00

LOJ (V<=1200, E<=120000):
Mine: 56 points
Mine w/ scaling: 88 points
Emaxx dinic: 44 points
LOJ fastest: AC, 1176 ms total, #1 on leaderboard

Synthetic Test(N=2000, E=4e6):
Mine: 66.56
Emaxx dinic: 247.40
Mine w/ scaling: 0.14
LOJ fastest: 0.10

Number of non-whitespace characters:
Mine: 1164
Emaxx dinic: 1326
LOJ fastest: 2228

## Conclusions:

First off, HLPP seems to be superior to Dinic's in terms of runtime, so perhaps the real lesson is that people should be using HLPP. Also, I hope to provide a solid well-tested implementation of Dinic's with scaling that people can use. Other than that, however, I hope it's clear that Dinic's with scaling really does have a significant improvement in runtime (approximately a factor of V), and that translates to significantly better performance on real problems.

PS: If anyone has a flow implementation they think is better than mine, I'd like to see it :^)

• +149

 » 11 months ago, # |   0 Auto comment: topic has been updated by Chilli (previous revision, new revision, compare).
 » 11 months ago, # | ← Rev. 4 →   +29 I've just been made aware of dacin21'simplementation, and although I haven't had the time to properly compare it yet — it does seem to be fairly fast — 0.05 on SPOJ, 0.03 on VNSPOJ, and 2000 ms on LOJ. It also seems to be reasonably short — seems very impressive.
•  » » 11 months ago, # ^ |   +21 It's not really my own implementation, I think I got it from the kactl codebook. If I recall correctly, the only thing I changed in this version was fixing a runtime error when the source node is isolated. (I also have a modified version of it for fast bipartite matching.)I'll try some of my other codes on the LOJ problem and see how they fare. (I couldn't really benchmark on the SPOJ problem anymore, as the runtime were really low. Moving to static memory allocation reduces the runtime from $0.05$ to $0.02$ for example. I didn't want to do that in my codebook for obvious reasons.)
•  » » » 11 months ago, # ^ |   +13 Kactl notebook truly is a gold mine of great implementations...
•  » » » 10 months ago, # ^ |   +11 Does push-relabel provide a noticeable speed improvement in bipartite matching over Hopcroft-Karp?
•  » » » » 10 months ago, # ^ |   +13 From my experience (from a few problems I solved) it does, although I have never benchmarked it (I never wrote Hopcroft-Karp actually). My $0.00$ on https://www.spoj.com/problems/MATCHING/ uses Push-Relabel.One thing to note for matching problem is that you should pick the vertex with the lowest label to discharge (for flow you pick the one with the highest label) and do global relabeling, both to get the theoretical $O(\sqrt{V} E)$ bound and for it to run fast in practice. There were also quite a few constant-factor optimizations, such as keeping all the excess on one side of the graph.
•  » » » » » 10 months ago, # ^ | ← Rev. 3 →   0 Do you have a source for the fact that "lowest global relabel" gives $O(\sqrt{V}E)$? I couldn't find anything from the google search I just did.Just found it on a comment from Dacin21: https://codeforces.com/blog/entry/58048?#comment-417654
 » 11 months ago, # | ← Rev. 2 →   +4 nice
 » 11 months ago, # |   +24 The LOJ submissions was really interesting to have a look at. Here's my quick analysis of it: It uses the highest label push-relabel with buckets. Edges are stored in a 2D adjacency list. (The alternative is to compress the graph into a 1D list.) It uses gap relabeling. Gap relabeling is optimized by maintaining all nodes (active or inactive) in an additional bucket datastructure. This is crucial on SPOJ, where repeated gab-relabellings are the bottleneck. (Normal gap relabeling runs in $O(n)$ every time it occurs, with this is O(n) in total. This does not matter on LOJ that much, as the graphs there are a lot denser.) This optimization was new to me. It only computes the min-cut and not the full max flow, as vertices with heights $\geq n$ are ignored. This only works on tasks where you only need the total flow and not the flow on specific edges. It does not maintain the current edge if each node. This leads to more arc-scans on LOJ. (The advantage of this is that you don't need a second pass to relabel, but the increase in arc-scans is also significant.) It does a single global relabeling step at the beginning. This helps in many cases. If we connect every node to both the source and the sink with edges of capacity $1$, then this won't help at all. With the improved gap relabeling, I now also got my $0.00$ on SPOJ, but getting that is rather inconsistent (I submitted a few times and the runtime varied between $0.00$, $0.01$ and $0.02$.)
•  » » 11 months ago, # ^ |   0 Thanks for the analysis of that implementation! I've noticed that for my implementation of Dinic's, storing the edges in a 2d adjacency list performs better than in a flat structure.The single global relabeling is interesting too — I've always thought it was something you do repeatedly. Is there any way to test its performance without that single global relabeling at the beginning? Removing it for me causes the code to fail.
•  » » » 11 months ago, # ^ |   +5 Yeah, running global updates every $\sim 4 n$ relabels is what the papers suggest. With some heuristics in the global relabeling, I got the following runtimes on LOJ (SPOJ is still $0.00$) (with static memory)$\require{siunitx}$ $\begin{array}[lr]\\ \text{global updates} & \text{runtime}\\ \hline \text{yes} & 947 \,\mathrm{ms}\\ \text{initial only} & 1334\,\mathrm{ms}\\ \text{no} & 1465\,\mathrm{ms} \end{array}$so it does indeed help. I also tried partial augmentation, but that was a lot slower. (Maybe I understood it wrong, the notes about it "The Partial Augment-Relabel Algorithm for the Maximum Flow Problem" are a bit vague about the multi-push operation.)
•  » » » » 11 months ago, # ^ | ← Rev. 2 →   0 Well, I get slightly less drastic results, but it does seem to be superior to do global updates every $4n$ relabels. One note, I'm counting each gap relabel as multiple relabels, while you seem to count them as just one — counting them as multiple relabels seems to help for me.With global updates, I get a best of $913$ ms, and with only an initial update, I get a best of $1050$ ms.
 » 11 months ago, # | ← Rev. 2 →   +16 If anyone wants to use the really fast implementation from LOJ in their library, I simplified it and cleaned up some stuff here: https://gist.github.com/Chillee/ad2110fc17af453fb6fc3357a78cfd28#file-hlpp-cppIt's currently the fastest/tied for 0.00 seconds on all 3 of SPOJ FASTFLOW, VNSPOJ FASTFLOW, and LOJ flow.It comes in at 87 lines and 1687 non whitespace characters, so decently long, but about 30% shorter than the original implementation and a lot easier to understand for me.
•  » » 11 months ago, # ^ | ← Rev. 2 →   0 There is an error in the code. In line 41, you try to get the back element from an empty queue.
•  » » » 11 months ago, # ^ |   0 Well... Perhaps undefined behavior is the correct word. At least in my compiler (gcc c++14), it returns a reference to the last element that was in the queue.I do agree that it's pretty jank(and should be changed), but it's not an "error".
•  » » » » 11 months ago, # ^ |   0 Sorry, I meant undefined behavior too. It surprisingly doesn't work on my compiler and gives a runtime error.
•  » » » » » 11 months ago, # ^ |   +11 Alright, I fixed it. Thanks for pointing out that it fails.
 » 10 months ago, # |   0 Hi, thanks for the nice implementation first! In SPOJ, I use my java version based on your hlpp code, and I got a rte(nzec). And when i use 2 * MAXN to create these arrays and vectors then I got an ac. So can you show your spoj submission? I want to check whether there are errors in my java implementation.
•  » » 10 months ago, # ^ |   0
•  » » » 10 months ago, # ^ |   0 Seems maxV + 5 will solve this issue. Thanks for sharing.


text:
## LuxRender official feature request thread

General Project and community related discussion.

Moderator: coordinators

## attention: please vote only one time, vote are registered as image

physical sky clouds texturing TOT. VOTES = 42 viewtopic.php?f=13&t=4373#p45465
19
23%
exif data TOT. VOTES = 10 viewtopic.php?f=13&t=4373&start=10#p45498
8
10%
camera color correction TOT. VOTES = 15 viewtopic.php?f=13&t=4373&start=10#p45544
6
7%
Face rendering behavior TOT. VOTES = 8 viewtopic.php?f=13&t=4373&start=10#p45554
3
4%
Media precedence TOT. VOTES = 13 viewtopic.php?f=13&t=4373&p=45779#p45777
7
8%
Multiple UVs per mesh TOT. VOTES = 27 viewtopic.php?f=13&t=4373&p=45779#p45777
14
17%
full motion blur TOT. VOTES = 28 viewtopic.php?f=13&t=4373&start=30#p45895
13
16%
Ptex Per-Face Texture Map. TOT. VOTES = 21 viewtopic.php?f=13&t=4373&start=30#p45981
13
16%

### Re: LuxRender official feature request thread

jeanphi wrote:The light group processing is done before the imaging pipeline (which includes gamma correction) this would be a rather intrusive change, and I'm not sure it'd bring much in practice. If you want your sky to be more blueish, why don't you apply a blue filter on it?

I think that we max user can't add a filter to the sky, LuxMax is different from LuxBlend.... we have control only about gain, sun size and turbidity.
probably to add a filter to the sky we need it to have as shader... as shader we can edit it in material editor... but max shaders aren't compatible with luxrender. I will pinpoint to Hedphelym to add the missing value A, B, C, D and E parameters as "float aconst", "float bconst" etc. that's way I asked for this small tool... a small tool to find as sub tool like the color piker for each created light layer.
I'm thinking about to ask to Hedphelym if there is a way to create an scripted external material and filter editor for luxrender.... similar to the one used from luxblend user.

jeanphi wrote:
patro wrote:a small preview render window in light group to see the result of all added and edited effects... before to waste time by adding it to the whole render.

This I really don't get, light group feedback is interactive in the whole render.

Jeanphi

I had think that this was impossible without an intrusive change... but lux have the update button in render UI. other render have this small window... and have also the update button. anyway it could be an idea for the future!

patro

Posts: 1798
Joined: Fri Feb 29, 2008 9:06 pm
Location: mount Etna

### Re: LuxRender official feature request thread

patro wrote:I think that we max user can't add a filter to the sky, LuxMax is different from LuxBlend....

I was just hinting to use the sky lightgroup color filter.

Jeanphi
jeanphi

Posts: 6577
Joined: Mon Jan 14, 2008 7:21 am

### Re: LuxRender official feature request thread

jeanphi wrote:
patro wrote:I think that we max user can't add a filter to the sky, LuxMax is different from LuxBlend....

I was just hinting to use the sky lightgroup color filter.

Jeanphi

ah OK. like I do here: viewtopic.php?f=36&t=4595&start=10#p48135
I understand now... how to use correctly the color piker in light group. thanks Jeanphi to let me think about....
see attached image as example
normal sky

wrong cruise poition... the horizon is near by fully override

right cruise position... the horizon is near by maintained.

but a contrast or hue-saturation tool will increase the blue from the sky(generic color)... by increasing difference between the sky and the transition part like the horizon by don't adding more color to the horizon/shadowed zone.

patro

Posts: 1798
Joined: Fri Feb 29, 2008 9:06 pm
Location: mount Etna

### Re: LuxRender official feature request thread

hi Jeanphi,

please could you add as soon as possible the "environment background mapping function" - tilling, rotating etc. and also to avoid that the background image will loaded flipped?
thank you

Patro

patro

Posts: 1798
Joined: Fri Feb 29, 2008 9:06 pm
Location: mount Etna

### Re: LuxRender official feature request thread

patro wrote:please could you add as soon as possible the "environment background mapping function" - tilling, rotating etc. and also to avoid that the background image will loaded flipped?

Well, as soon as possible might not be soon I currently have very few time to spend developing LuxRender. I'd like to get out v0.7.1 first otherwise it will never be released.

Jeanphi
jeanphi

Posts: 6577
Joined: Mon Jan 14, 2008 7:21 am

### Re: LuxRender official feature request thread

patro wrote:please could you add as soon as possible the "environment background mapping function" - tilling, rotating etc. and also to avoid that the background image will loaded flipped?

I'm not exactly sure what you're asking for here, could you please elaborate? Infinite light (used for HDR environment maps) already supports rotation, and tiling makes little sense to me.
May contain traces of nuts.

Lord Crc

Posts: 4454
Joined: Sat Nov 17, 2007 2:10 pm

### Re: LuxRender official feature request thread

Lord Crc wrote:I'm not exactly sure what you're asking for here, could you please elaborate? Infinite light (used for HDR environment maps) already supports rotation, and tiling makes little sense to me.

as I know isn't possible to rotate an image used as environment... around Z axis.
you can set an exr an jpg.. but I can't see a spinner to rotate the image like I want.
and I asked more than once about... and I got always the same answer: not supported at time!
I'm talking about a latlong image... to use as environment illumination and as background... if I load in 3dsmax viewport an latlong image, I can rotate the image... so to display the portion that I want but when I render the scene... the image isn't rotated like in my viewport and the image is rendered also flipped.
you are saying me that rotation is already supported.. OK... below is the exported code for the background environment image..please could you add the missing code part for rotate it?

# Enviromentmap goes here
AttributeBegin
LightGroup "environment"
LightSource "infinite"
"string mapping" ["latlong"]
"float gamma" [1.0]
"float gain" [1.0]
AttributeEnd

thanks

PS: about tailing... consider the use of a small checker image...

patro

Posts: 1798
Joined: Fri Feb 29, 2008 9:06 pm
Location: mount Etna

### Re: LuxRender official feature request thread

patro wrote:you are saying me that rotation is already supported.. OK... below is the exported code for the background environment image..please could you add the missing code part for rotate it?

Code: Select all
#Env userprop F:\CAD\mappe\hrdi dds\StudioHDR\small 03.exr# Enviromentmap goes hereAttributeBegin # rotate image 180 degrees around Z axisRotate 180.0  0 0 1 LightGroup "environment" LightSource "infinite" "string mapping" ["latlong"] "string mapname" ["F:\\CAD\\mappe\\hrdi dds\\StudioHDR\\small 03.exr"] "float gamma" [1.0] "float gain" [1.0] AttributeEnd

This is already exposed in LuxBlend. Rotations are supported, I tried flipping around the standard axes by scaling by -1, which seems to work, however I don't know the exact limits for which kind of transformations are allowed. Perhaps jeanphi can elaborate?
May contain traces of nuts.

Lord Crc

Posts: 4454
Joined: Sat Nov 17, 2007 2:10 pm

### Re: LuxRender official feature request thread

I can't understand! I asked more than once for something like this:

# rotate image 180 degrees around Z axis
Rotate 180.0 0 0 1

but the answer wasn't: rotation supported and tailing no.

thanks a lot... I will pinpoint this to Hedphelym to add this part to the exporter.

patro

Posts: 1798
Joined: Fri Feb 29, 2008 9:06 pm
Location: mount Etna

### Re: LuxRender official feature request thread

patro wrote:I can't understand! I asked more than once for something like this

Hmm only I can think of then is that you must have been misunderstood. At least we got it sorted now

Please note that the rotation/transformation is applied independently of the mapping (ie angular, latlong).

edit: that said, tiling is not supported. Again I can't see why this is useful, it would give incorrect/weird lighting.
May contain traces of nuts.

Lord Crc

Posts: 4454
Joined: Sat Nov 17, 2007 2:10 pm

PreviousNext

format: plain


text:
1. Calculation functions
2. All functions
3. Financial functions
4. DURATION

You can use the DURATION function to calculate the Macauley duration for an assumed parity value of 100 monetary units.

The Macauley duration is the weighted average maturity of cash flows. That is, the weighted average distance to payment. It's used to measure a bond price's response to changes in yield. A higher Macauley duration value indicates a riskier investment.

## Syntax

DURATION(Settlement, Maturity, Rate, Yield, Frequency [, Basis])

## Arguments

 Argument Data type Description Settlement (required) Date The bond settlement date: The date the bond is traded to the buyer. Maturity (required) Date The bond maturity date: The date when the bond expires. Rate (required) Number The bond annual coupon date. Yield (required) Number The bond annual yield. Frequency (required) Number The number of coupon payments per year. Enter: 1 for annual 2 for semi-annual 4 for quarterly Basis Number The basis determines how many days exist in a year. A full year has: 360 days when basis US (NASD) 30/360, Actual/360, and EUR 30/360 are used 365 days when basis Actual/365 is used 365 or 366 days when Actual/Actual is used US 30/360 is the default basis for COUPDAYS. It can also be specified by entering 0. To use a different type of day count basis, enter: 1 for Actual/Actual 2 for Actual/360 3 for Actual/365 4 for European 30/360 Learn about the conventions used to calculate the day count for basis.

The DURATION function returns a number.

The Macauley duration is calculated with the following formula:

$Duration = \frac { \sum_{t=1}^{T} \frac {tC} {(1+y)^{t}} + \frac{TF}{(1+t)^{T}} } { P }$

Where:

• C is coupon
• y is yield
• F is face value
• P is price, inclusive of accrued interest
• T is number of periods

## Constraints

• The settlement and maturity dates must be valid dates between 01/01/1900 and 12/31/2399.
• The maturity date must be later than the settlement date.
• The rate and yield must be positive or zero.
• The frequency must be either 1 (annual), 2 (semi-annual), or 4 (quarterly).
• The basis, when specified, must be either 0 (US 30/360), 1 (Actual/Actual), 2 (Actual/360), 3 (Actual/365), or 4 (EUR 30/360).

DURATION

## Examples

This example shows a Macauley duration calculation that specifies a basis.

 Formula Description Result DURATION(DATE(2018, 1, 15), DATE(2021, 1, 15), 0.12, 0.1, 1, 4) The example has a: settlement date of 01/15/2015 maturity date of 01/15/2018 rate of 0.12 (12%) yield of 0.1 (10%) frequency of 1 (annual) basis of 4 (European 30/360) 2.6976811

This example shows a Macauley duration calculation that does not specify a basis. As a result, the basis defaults to US 30/360.

 Formula Description Result DURATION(DATE(2018, 1, 15), DATE(2021, 1, 15), 0.12, 0.1, 4) The example has a: settlement date of 01/15/2018 maturity date of 01/15/2021 rate of 0.12 (12%) yield of 0.1 (10%) frequency of 4 (quarterly) 2.5760086

Disclaimer

We may update our documentation occasionally, but will only do so in a way that does not negatively affect the features and functionality of the Anaplan service.


text:
Goldbach’s Comet with Numba and Datashader

This Python post is about computing and plotting Goldbach function. It is using some very basic mathematical knowledge, nothing fancy! The point is to do some computations with Numba. We also use Datashader to perform some efficient plotting.

Here is the definition of the Goldbach function from wikipedia:

The function $g ( E )$ is defined for all even integers $E > 2$ to be the number of different ways in which E can be expressed as the sum of two primes. For example, $g ( 22 ) = 3$ since 22 can be expressed as the sum of two primes in three different ways ( 22 = 11 + 11 = 5 + 17 = 3 + 19).

Note that for Goldbach's conjecture to be false, there must be $g(E) = 0$ somewhere 🤯.

Anyway, here are the steps used in this post to compute Goldbach function:

• Define a maximum positive integer $n$.
• For each natural number smaller or equal to $n$, build a quick way to check if it is a prime or not. In order to do that, we are going to create a boolean vector using the sieve of Eratosthenes.
is_prime_vec = generate_is_prime_vector(n)
• for each even number $E$ smaller or equal to $n$, compute $g(E)$ by counting the number of cases where $E-p$ is prime for all primes $p$ not larger than $E/2$ [ [$p$, $E-p$] is a Goldbach partition of $E$ if they are both primes].
g_vec = compute_g_vector(is_prime_vec)

Imports

import numpy as np
import pandas as pd
import primesieve
from numba import jit, njit, prange
import matplotlib.pyplot as plt
import perfplot
from datashader import transfer_functions as tf
from colorcet import palette

Package versions:

    Python    : 3.9.6
pandas    : 1.3.0
primesieve: 2.3.0
matplotlib: 3.4.2
numpy     : 1.21.1
colorcet  : 2.0.6
perfplot  : 0.9.6
numba     : 0.53.1

Sieve of Eratosthenes with Numba

This code is inspired from wikipedia and geeksforgeeks. It flags all primes not greater than $n$. The algorithm includes a common optimization, which is to start enumerating the multiples of each prime $p$ from $p^2$.

@jit(nopython=True)
def generate_is_prime_vector(n: int) -> np.ndarray:
"""Sieve of Eratosthenes.

Parameters
----------
n : int
n is the largest number that we want to flag ad prime or not.

Returns
-------
is_prime_vec : ndarray
an array of type bool (i is prime if is_prime_vec[i] is True).
"""

# initialize all entries as True, except 0 and 1
is_prime_vec = np.ones(n + 1, dtype=np.bool_)
is_prime_vec[0:2] = False

# loop on prime numbers
p = 2
while p * p <= n:
if is_prime_vec[p] == True:
# Update all multiples of p, starting from p * p
for i in range(p * p, n + 1, p):
is_prime_vec[i] = False
p += 1

return is_prime_vec

Let's run it with a small value of n:

n = 11
generate_is_prime_vector(n)
array([False, False,  True,  True, False,  True, False,  True, False,
False, False,  True])

We can get the list of primes from the is_prime vector using np.nonzero:

np.nonzero(generate_is_prime_vector(n))[0]
array([ 2,  3,  5,  7, 11])

The list of primes not greater then n can also be computed using the primesieve library:

primesieve.primes(n)
array('Q', [2, 3, 5, 7, 11])

Let's check that the list of primes found is correct with a larger value of n:

%%time
n = 1_000_000
is_prime_vec = generate_is_prime_vector(n)
prime_list_1 = list(np.nonzero(is_prime_vec)[0])
prime_list_2 = primesieve.primes(n)
np.testing.assert_array_equal(prime_list_1, prime_list_2)
CPU times: user 37.1 ms, sys: 69 µs, total: 37.2 ms
Wall time: 36.5 ms

Elapsed time

This generate_is_prime_vector function is rather efficient compared to the next step of computing $g(E)$ for all $E$s. This classical implementation of the sieve of Eratosthenes is supposed to be $O(n log log n)$.

out = perfplot.bench(
setup=lambda n: n,
kernels=[
lambda n: generate_is_prime_vector(n),
],
labels=["generate_is_prime_vector(n)"],
n_range=[10 ** k for k in range(1, 10)],
)
out
┏━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓
┃ n          ┃ generate_is_prime_vector(n) ┃
┡━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩
│ 10         │ 6.3e-07                     │
│ 100        │ 7.070000000000001e-07       │
│ 1000       │ 1.7830000000000001e-06      │
│ 10000      │ 1.2794000000000002e-05      │
│ 100000     │ 0.000149319                 │
│ 1000000    │ 0.002291595                 │
│ 10000000   │ 0.081507009                 │
│ 100000000  │ 1.119336045                 │
│ 1000000000 │ 11.769184520000001          │
└────────────┴─────────────────────────────┘

ms = 10
fig = plt.figure(figsize=(10, 10))
plt.loglog(
out.n_range,
1.0e-9 * out.n_range * np.log(np.log(out.n_range)),
"-",
label="$O(n \; log \; log \; n)$",
)
plt.loglog(
out.n_range, out.timings_s[0], "v-", ms=ms, label="generate_is_prime_vector(n)"
)
plt.legend()
plt.grid(True)
_ = ax.set_ylabel("Runtime [s]")
_ = ax.set_xlabel("n")
_ = ax.set_title("Timings of generate_is_prime_vector")

The timings of the same function in pure Python is not shown. It is usually slower by a factor 10 or 100 [100 is more likely than 10 for this kind of CPU bound function].

Find the number of prime pairs for a given even number E

Now we show how $g(E)$ can be computed for a given value of $E$:

# E is the even number that we want to decompose into prime pairs
E = 22
assert E % 2 == 0

E_half = int(0.5 * E)

# we generate the boolean prime vector for all the integer up E
is_prime_vec = generate_is_prime_vector(E)

# initialization
i = 2  # loop counter (we know that 0 and 1 are not primes)
count = 0  # number of prime pairs

# we loop over all the prime numbers smaller than or equal to half of num
while i <= E_half:
if is_prime_vec[i] and is_prime_vec[E - i]:
print(f"({i:2d}, {E-i:2d})")
count += 1
i += 1
print(f"{count} prime pairs")
( 3, 19)
( 5, 17)
(11, 11)
3 prime pairs

Loop over all even numbers E not larger than n

Now we apply the previous process to all $E$s smaller or equal to $n$. We only compute is_prime_vec once and use it for all the evaluations of $g(E)$. Note that in the following compute_g_vector function, the outer loop has a constant step size of 1, in order to later use Numba prange, which only supports this unit step size. This means that we loop on contiguous integer values of $E/2$ instead of even values of $E$.

@jit(nopython=True)
def compute_g_vector(is_prime_vec: np.ndarray) -> np.ndarray:
"""Evaluate the Goldbach function.

This evatuates the Goldbach function for all even integers E not
larger than the largest integer from the argument is_prime.

Parameters
----------
is_prime_vec : ndarray
an array of type bool (i is prime iif is_prime[i] is True).

Returns
-------
g_vec : ndarray
an array of type uint (g_vec[i] corresponds to g(E) where E = 2i).
"""

n_max = len(is_prime_vec) - 1
n_max_half = int(0.5 * n_max) + 1

g_vec = np.empty(n_max_half, dtype=np.uint)

for i in range(n_max_half):
count = 0
E = 2 * i
for j in range(2, i + 1):
if is_prime_vec[j] and is_prime_vec[E - j]:
count += 1
g_vec[i] = count

return g_vec

The $i$-th value of g_vec correponds to $g(2 i)$ with $i \geq 0$:

i E g_vec[i]
0 0 0
1 2 0
2 4 1
3 6 1
4 8 1
5 10 2

and so on... We can check the values of $g$ at least for some for some small values of $E$:

n = 56
is_prime_vec = generate_is_prime_vector(n)
g_vec = compute_g_vector(is_prime_vec)
g_vec_ref = [0, 0, 1, 1, 1, 2, 1, 2, 2, 2, 2, 3, 3, 3, 2, 3, 2, 4, 4, 2, 3, 4,
3, 4, 5, 4, 3, 5, 3]
np.testing.assert_array_equal(g_vec, g_vec_ref)

First plot of the comet

%%time
n = 10_000
is_prime_vec = generate_is_prime_vector(n)
g_vec = compute_g_vector(is_prime_vec)
g_df = pd.DataFrame(data={"g": g_vec}, index=2 * np.arange(len(g_vec)))
g_df = g_df[g_df.index > 2]  # The function g(E) is defined for all even integers E>2
g_df.head(2)
CPU times: user 59.9 ms, sys: 150 µs, total: 60.1 ms
Wall time: 60.5 ms

g
4 1
6 1
ax = g_df.plot(style=".", ms=5, alpha=0.5, legend=False, figsize=FS)
_ = ax.set(
title="Goldbach's comet",
xlabel="E",
ylabel="g(E)",
)
ax.autoscale(enable=True, axis="x", tight=True)
ax.grid(True)

Parallel loop

We are now going to parallelize compute_g_vector just by using the Numba njit decorator and prange on the outer loop. The vector g_vec is shared across threads but each thread is only writing into its own partition.

@njit(parallel=True)
def compute_g_vector_par(is_prime_vec: np.ndarray) -> np.ndarray:
"""Evaluate the Goldbach function.

This evatuates the Goldbach function for all even integers E not
larger than the largest integer from the argument is_prime.

Parameters
----------
is_prime_vec : ndarray
an array of type bool (i is prime iif is_prime[i] is True).

Returns
-------
g_vec : ndarray
an array of type uint (g_vec[i] corresponds to g(E) where E = 2i).
"""

n_max = len(is_prime_vec) - 1
n_max_half = int(0.5 * n_max) + 1

g_vec = np.empty(n_max_half, dtype=np.uint)

for i in prange(n_max_half):
count = 0
E = 2 * i
for j in range(2, i + 1):
if is_prime_vec[j] and is_prime_vec[E - j]:
count += 1
g_vec[i] = count

return g_vec

Here we check that the sequential and parallel versions return the same values:

g_vec_par = compute_g_vector_par(is_prime_vec)
np.testing.assert_array_equal(g_vec, g_vec_par)

Elapsed time

Now we compare the running time of the sequential and parallel versions of compute_g_vector. The code is running on a laptop with 8 cores [Intel i7-7700HQ CPU @ 2.80GHz].

out = perfplot.bench(
setup=lambda n: generate_is_prime_vector(n),
kernels=[
lambda is_prime_vec: compute_g_vector(is_prime_vec),
lambda is_prime_vec: compute_g_vector_par(is_prime_vec),
],
labels=["jit", "njit"],
n_range=[10 ** k for k in range(1, 6)],
)
out
┏━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━┓
┃ n      ┃ jit                   ┃ njit                   ┃
┡━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━┩
│ 10     │ 8.480000000000001e-07 │ 6.506e-06              │
│ 100    │ 2.271e-06             │ 7.0370000000000006e-06 │
│ 1000   │ 0.000224464           │ 6.1258e-05             │
│ 10000  │ 0.033075403           │ 0.010347824            │
│ 100000 │ 2.923607557           │ 1.1248999480000001     │
└────────┴───────────────────────┴────────────────────────┘

Basically, running the parallel version on 8 cores divides the running time by 3.

ms = 10
fig = plt.figure(figsize=(10, 10))
plt.loglog(
out.n_range,
1.0e-9 * np.power(out.n_range, 2),
"-",
label="$O(N^2)$",
)
plt.loglog(out.n_range, out.timings_s[0], "v-", ms=ms, label="jit")
plt.loglog(out.n_range, out.timings_s[1], "^-", ms=ms, label="njit")
plt.legend()
plt.grid("on")
_ = ax.set_ylabel("Runtime [s]")
_ = ax.set_xlabel("N")
_ = ax.set_title("Timings of compute_g_vector")

If we run the parallel version with $n=1e6$, we might estimate the running time of $n=1e7$ and see if this is affordable.

%%time
n = 1_000_000
is_prime_vec = generate_is_prime_vector(n)
g_vec_par = compute_g_vector_par(is_prime_vec)
CPU times: user 6min 9s, sys: 428 ms, total: 6min 9s
Wall time: 1min 28s

So computing $g$ with $n = 1e7$ should take one or a few hours on the laptop executing compute_g_vector_par... Let's do that next.

Compute more data

%%time
n = 10_000_000
is_prime_vec = generate_is_prime_vector(n)
CPU times: user 162 ms, sys: 0 ns, total: 162 ms
Wall time: 68.7 ms

%%time
g_vec_par = compute_g_vector_par(is_prime_vec)
g_df = pd.DataFrame(data={"E": 2 * np.arange(len(g_vec_par)), "g": g_vec_par})
g_df = g_df[g_df.E > 2]  # The function g(E) is defined for all even integers E>2
CPU times: user 8h 48min 44s, sys: 558 ms, total: 8h 48min 44s
Wall time: 1h 49min 9s

A bit less than 2 hours...

cmap = palette["kbc"][::-1]
bg_col = "white"
height = 800
width = int(np.round(1.6180 * height))
%%time
cvs = ds.Canvas(plot_width=width, plot_height=height)
agg = cvs.points(g_df, "E", "g")
img = tf.set_background(img, bg_col)
img
CPU times: user 809 ms, sys: 23.9 ms, total: 833 ms
Wall time: 826 ms

We can clearly observe some dense lines in this "comet tail". In order to visualize this vertical distribution of prime pairs count, we are are going to normalize $g$. As explained on the wikipedia page:

An illuminating way of presenting the comet data is as a histogram. The function $g(E)$ can be normalized by dividing by the locally averaged value of g, gav, taken over perhaps 1000 neighboring values of the even number E. The histogram can then be accumulated over a range of up to about 10% either side of a central E.

%%time
g_df["g_av"] = g_df["g"].rolling(window=1000, center=True).mean()
g_df["g_norm"] = g_df["g"] / g_df["g_av"]
CPU times: user 170 ms, sys: 67.8 ms, total: 238 ms
Wall time: 230 ms

%%time
cvs = ds.Canvas(plot_width=width, plot_height=height)
agg = cvs.points(g_df, "E", "g_norm")
img = tf.set_background(img, bg_col)
img
CPU times: user 644 ms, sys: 27.9 ms, total: 672 ms
Wall time: 674 ms

We can now proceed to plot the histogram of the comet data, which will lead to some kind of cross section of the above plot:

ax = g_df["g_norm"].hist(bins=1000, alpha=0.5, figsize=(20, 10))
ax.grid(True)
_ = ax.set(
title="Histogram of the normalized Goldbach function",
xlabel="Normalized number of prime pairs",
ylabel="Number of occurrences",
)
_ = plt.xticks(np.arange(0.5, 3, 0.1))

Prime E/2 values only

Finally, we are going to isolate a part of the most dense line from the comet tail [for a normalized number of prime pairs around 0.66]. As explained in wikipedia:

Of particular interest is the peak formed by selecting only values of E/2 that are prime. [...] The peak is very close to a Gaussian form.

g_df["E_half"] = (0.5 * g_df["E"]).astype(int)
g_df_primes = g_df[g_df["E_half"].isin(np.nonzero(is_prime_vec)[0])]
%%time
cvs = ds.Canvas(plot_width=width, plot_height=height)
agg = cvs.points(g_df_primes, "E", "g_norm")
img = tf.set_background(img, bg_col)
img
CPU times: user 79 ms, sys: 0 ns, total: 79 ms
Wall time: 75.9 ms

ax = g_df_primes["g_norm"].hist(bins=500, alpha=0.5, figsize=(20, 10))
ax.grid(True)
_ = ax.set(
title="Histogram of the normalized Goldbach function for prime E/2 values only",
xlabel="Normalized number of prime pairs",
ylabel="Number of occurrences",
)


text:
# What can you do with this?

We interrupt our regularly scheduled long patch of radio silence to share with you an arresting, mathematically rich visual:

[source]

Thanks to Josh Kershenbaum for the tip.

UPDATE 7/17:

To clarify something: this image was originally created as an answer to a question, but I didn’t create a WCYDWT tag for that angle on it. To me, the image lands as (a) very beautiful; relatedly, (b) haunting, hard to get out of my head; therefore, (c) incipiently highly narrative – it is asking us to surround it with story, whether the original story conceived by the maker of the image or some other; and, lastly, (d) unavoidably mathematical. I don’t have a clear sense of the next move, but I do think this image opens a rich vein of something for a math teacher to use. The question “What can you do with this?” isn’t rhetorical: what’s your next move?

(See exchange with Dan Meyer below.)

# The Inconvenient Truth Behind Waiting for Superman and other stories

I just recently learned of an organization in NYC called the Grassroots Education Movement, which last Thursday premiered a documentary film with the awesome title The Inconvenient Truth Behind Waiting for Superman. They will apparently send you a copy for free; I just ordered mine.

Meanwhile, the city of New York continues to besiege its own public schools with budget cuts, looming layoffs, and a multi-year hiring freeze. (Having spent the year training 12 new teachers, let me not even get started on the hiring freeze.) Another thing that happened on Thursday was that East Side Community High School, a wonderful school on the Lower East Side where I used to teach and where the math teaching is strong enough that we placed four student teachers there this quarter, had its first fund-raiser. Like, big event, speakers, performances by students, pay to get in, as though it were a non-profit, carrying out its own civic mission and in need of private funding to do it, rather than a public school, charged with a civic mission by the state, which no longer sees fit to pay for it.

I missed both the documentary premiere and ESCHS’s fund-raiser because I was teaching the final class of a 3-session minicourse at Math for America on the fundamental theorem of arithmetic. Let me do a little reflecting on the execution:

At the end of the 2nd session, I gave participants about a half-hour to try to figure out something quite difficult. I attempted to scaffold this with some unobtrusive PCMI-style tricks in a previous problem set: sequences of problems with the same answer for a mathematically significant reason. It turned out not to be enough. There was high engagement the whole time, but no one seemed to be headed in my intended direction after that half-hour. On the other hand, that half-hour had made the group into a legit mathematical research community. What was afoot was a live process of trying things out, questioning, pressing on others’ logic, and generally behaving like research mathematicians. I was left with a dilemma. I had one session remaining. I wanted to protect that process, meaning I did not want to steal from them any of the deliciousness (or pain – also delicious) of the process they were in the middle of by offering them too much direction. But at the same time I felt I needed to guarantee that we would reach resolution. (Storytelling purposes.)

The solution I went with: I had them pick up in the final session where they left off, but I brought in a sequence of hints on little cut-up slips of paper. I tried to call them “idea-starters” as opposed to “hints” to emphasize that the game was you thinking on your own, and this is just to get you moving if you’re stuck, rather than I have a particular idea and I want you to figure out what it is, but I don’t think I was consistent with this, and I think they pretty much all called them “hints,” and I don’t think it really mattered. They were in an order from least-obtrusive to most-directive. None of them were very directive. Most importantly, I told the participants that if they wanted to get one, they needed to decide this as a table. (There were 6 tables with 3-4 folks each.)

How this went: a) it preserved the sense of mathematical community. I do not think there was much of a cost to participant ownership of what they found out. b) People were actually pretty hesitant to use the “idea-starters.” Most of them went untouched. This would probably be different with a different audience. (High schoolers instead of teachers?) c) The “idea-starters” worked great, but very slowly. I planned to spend 45 min letting them work in this arrangement, but after 45 min, most of the groups were still deep in the middle of something. After over an hour, I asked two groups to present what they had, however incomplete, for the sake of a change of pace and the opportunity for cross-pollination of ideas between the tables. I had actually meant to do a lot more of this but had forgot to mention it at the beginning. I let everybody work for another 10-15 min while these groups laid out their presentations. By the time they presented, I realized that there wasn’t enough time left for everyone to really get back to work afterward, but in any case their ideas had gotten more fully developed in that 10 min. so they actually had pretty much figured out everything I had wanted them to. I presented the final link in the logical chain, just to fill in the picture, in the last 5 minutes. It was pretty satisfying to me to watch the presentations, except that it happened so late in the session. This for two reasons. One was that I would have ideally liked to have time to encourage the participants to interrogate the presenters more, but there wasn’t time for that. The other was that I had intended to spend the last half-hour with the participants consolidating their understanding of the argument by applying it to a new situation in which they didn’t know the outcome and it would tell them; but we didn’t have time for that either. I really feel a loss about that.

If I were to repeat it I think I would interrupt much earlier to have people present partial work. The cross-pollination of ideas might or might not accelerate the figuring-out process. Either way I think the change of pace would have been good for concentration. Also, I could have put some of the questions I used as “idea-starters” into the Session 2 problem sets, trying to move some of the combustion I got in session 3 into session 2. But these would both be experiments as well. I hope I get a chance to try them.

# Angle Sum Formulas: Request for Ideas

One of the student teachers I supervise is planning a lesson introducing the sine and cosine angle sum formulas. I wanted to give him some advice on how to make the lesson better – in particular, along the axes of motivation and justification – and realized that, never having taught precalculus, I barely had any! Especially re: justification. I basically understand these formulas as corollaries of the geometry of multiplication of complex numbers.[1] I have seen elementary proofs, but I remember them as feeling complicated and not that illuminating.

So: how do you teach the trig angle sum formulas? And in particular:

* How do you make them seem needed? (I offered my young acolyte the idea of asking the kids to find sin 30, sin 45, sin 60, sin 75 and sin 90 – with the intention of having them be slightly bothered by the fact that they can do all but sin 75.)

* Do you state the formulas or do you set something up to have the kids conjecture them? If the latter, how do you do it? How does it fly?

* How do you justify them? Do you do a rigorous derivation? Do you do something to make them seem intuitively reasonable? What do you do and how does it fly?

* Do you do them before or after complex numbers, and do you connect the two? If so, how do you do it and how does it fly?

Any thoughts would be much appreciated.

Thanks to John Abreu, who sent me the following in an email –

Please find attached a Word document with the proofs of the trig angle sum formulas. After opening the document you’ll see a sequence of 14 figures, the conclusions are obtained comparing the two of them in yellow. Also, I left the document in “crude” format so it’ll be easier for you to decide the format before posting.

I must say that the proofs/method is not mine, but I can’t remember where I learned them.

with an attachment containing the following figures (click to enlarge / for slideshow) –

As far as I can tell, the proof is valid for any pair of angles with an acute sum.

Notes:

[1]Let $z_1,z_2$ be two complex numbers on the unit circle, at angles $\theta_1,\theta_2$ from the positive real axis. Then $z_1=\cos{\theta_1}+\imath\sin{\theta_1}$ and $z_2=\cos{\theta_2}+\imath\sin{\theta_2}$, so by sheer algebra, $z_1z_2=(\cos{\theta_1}\cos{\theta_2}-\sin{\theta_1}\sin{\theta_2})+\imath(\cos{\theta_1}\sin{\theta_2}+\sin{\theta_1}\cos{\theta_2})$. On the other hand, the awesome thing about multiplication of complex numbers is that the angles add – the product $z_1z_2$ will be at an angle of $\theta_1+\theta_2$ from the positive real axis; thus it is equal to $\cos{(\theta_1+\theta_2)} + \imath\sin{(\theta_1+\theta_2)}$. This is QED for both formulas if you believe me about the awesome thing. Of course it usually gets proven the other way – first the trig formulas, then use this to prove angles add when you multiply. But I think of the fact about multiplication of complex numbers as more essential and fundamental, and the sum formulas as byproducts.

# Creating Balance III / Miscellany

The Creating Balance in an Unjust World conference is back! I went a year and a half ago and it was awesome. Math education and social justice, what more could you want?

If you’re in NYC and you’re around this weekend, it’s happening right now! I’m going to try to make it to Session 3 this afternoon. It’s at Long Island University, corner of Flatbush and DeKalb in Brooklyn, right off the DeKalb stop on the Q train. I heard from one of the organizers that you can show up and register at the conference. I’m not 100% sure how that works given that it’s already begun, but I am sure you can still go.

* * * * *

I’ve just had a very intense week.

I want to get some thoughts down. I’m going to try very hard to resist my natural inclinations to a) try to work them into an overall narrative, and b) take forever doing it. Let’s see how I do.

(Ed. note: apparently not very well.)

I. Last spring I wrote

20*20 is 400; how does taking away 2 from one of the factors and 3 from the other affect the product? We get kids thinking hard about this and it would support the most contrivance-free explanation for why (neg)(neg)=(pos) that I have ever seen.

Without going into contextual details, let me just say that if you try to use this to actually develop the multiplication rules in a 1-hour lesson, all that will happen is that you will be dragging kids through the biggest, clunkiest, hardest-to-swallow, easiest-to-lose-the-forest-for-the-trees, totally-mathematically-correct-but-come-now model for signed number multiplication that you have ever seen (and this includes the hot and cold cubes). This idea makes sense for building intuition about signed numbers slowly, before they’re an actual object of study. It does not make any sense at all for teaching a one-off lesson explicitly about them. (Yes, the hard way. I totally knew this five months ago – what was I thinking?)

II. I gave a workshop Wednesday night, for about 35 experienced teachers, entitled “Why Linear Algebra Is Awesome.” The idea was to reinterpret the Fibonacci recurrence as a linear transformation and use linear algebra to get a closed form for the Fibonacci numbers. Again, without going into details –

I gave a problem set to make participants notice that the transformation we were working with was linear. I used those PCMI-style tricks like giving two problems in a row that have the same answer for a mathematically significant reason. This worked totally well. Here is the problem set:

Oops I guess I failed to avoid going into details. Anyway, the question was about how to follow this up. I went over 1-4 with everyone (actually, I had individual participants come up to the front for #3 and 4) at which point the only thing I really needed out of this – the linearity of the transformation – had been noticed by pretty much the whole room. One participant had gotten to #9 where you prove it, and I had her go over her proof.

I think this was valueless for the group as a whole. The proof was just a straight computation. You kind of have to do it yourself to feel it at all. It was such a striking difference watching people work on the problem set and have all these lightbulbs go off, vs. listening to somebody prove the thing they’d noticed. It almost seemed like people didn’t see the connection between what they’d noticed and what just got proved. I told them to take 5 minutes and discuss this connection with their table, but I got the feeling that this instruction was actually further disorienting for some participants.

I’m trying to put the experience into language so I get the lesson from it.

It’s like, there was something uninspired and disconnected about watching somebody formally prove the result, and then afterward trying to find the connection between the proof and the observation. Now that I write this down, clearly that was backward. If I wanted the proof (which was really just a boring calculation) to mean anything, especially if I wanted it to be at all engaging to watch somebody else do the proof, we needed to be in suspense about whether the result was true; either because we legitimately weren’t sure, or because we were pretty sure but a lot was riding on it.

This is adding up to: next time I do it, feel no need to prove the linearity. Let them observe it from the problem set and articulate it, but if there is no sense of uncertainty about it, this is enough. Later in the workshop, when we use it to derive a closed form for the Fibonacci numbers, now a lot is riding on it. If it feels right, we could take that moment to make sure it’s true.

III. As I work on my teacher class, something that’s impressing itself upon me for the first time is that definitions are just as important as proofs. What I mean by this is two things:

a) It makes sense to put a real lot of thought into motivating a course’s key definitions,

and maybe even more importantly,

b) Students of math need practice in creating definitions. You know I think that creating proofs is an underdeveloped skill for most students of math; it strikes me that creating definitions might be even more underdeveloped.

Definitions are one of the most overtly creative products of mathematical work, but they also solve problems. Not in quite the same sense that theorems do – they don’t answer precisely stated questions. But they answer an important question nonetheless – what do we really mean? And to really test a definition, you have to try to prove theorems with it. If it helps you prove theorems, and if the picture that emerges when you prove them matches the image you had when you started trying to make the definition, then it is a “good” definition. (This got clear for me by reading Stephen Maurer’s totally entertaining 1980 article The King Chicken Theorems.)

Anyway this adds up to an activity to put students through that I’ve never explicitly thought about before, but now find myself building up to with my teacher class:

a) Pose a definitional problem. Do a lot of work to make the class understand that we have an important idea at hand for which we lack a good definition.

b) Make them try to create a definition.

c) If they come up with something at all workable, have them try to use it to prove something they already believe true. I’ve often talked in the past about how trying to prove something you already believe true is very difficult, and that will be a problem here. However, unlike in the cases I had in mind (e.g. a typical Geometry “proof exercise”), this situation has the necessary element of suspense: does our definition work?

If they don’t come up with something workable, maybe give them a not entirely precise definition to try out.

d) Refine the definition based on the experience trying to use it to prove something.

I’ll let you know how it goes. I’m excited about it because it mirrors the process that advances mathematics as a discipline. But I expect to have a much better sense of its usefulness once I’ve given it an honest whirl.

Two things I forgot to say last night.

I. The reason I’m excited about the idea of having my class use its own self-made definitions to try to prove things is not just, or even primarily, because it will help them realize the inadequacies in their definitions. Although it will do that for sure. Even more than that, it seems to me the perfect way to support them in coming up with better definitions. This is what happened to Cauchy: he defined the limit verbally and a little vaguely, but then when he actually tried to use his definition to prove things, he started writing down precise inequalities. He didn’t have a teacher around to point out that this meant he should probably revise his definition, but my class does.

II. Yesterday when I asked my class to try to make a precise definition for what it means to converge, or for something to have a limit, some of them who took real analysis long ago began accessing this knowledge in an incomplete way. They started to talk about $\epsilon$ and $\delta$, but in vague, uncertain terms. It looked as though others might possibly accept the half-remembered vagueries because they seemed like they might be the “this is supposed to be the answer” answer. I had to prevent this. (The danger would have been even greater if these participants had correctly and confidently remembered the definition.) I stepped in to the conversation to say, yes, that thing you’re half-remembering is my objective, but what’s going to make you understand it so you never forget it again is to fight till you’re satisfied we’ve captured the meaning of convergence. You can either fight with the definition you half-remember or you can fight to build a new definition, but you have to go through your dissatisfaction to get there. You have to air all this dissatisfaction.

Afterward, I thought of a better language. I’ll give this to them next time.

Dissatisfaction is the engine that created analysis. This content, more than any other content, is both confusing and pointless if you bury your dissatisfaction rather than allowing it to thrive and be answered. The primary virtue of the tools of analysis is that they are satisfying. Only if you bring forth your dissatisfaction will this content have a chance to show you its value. So. Honor your dissatisfaction. It is the engine that will move us forward.

# Over the Course of an Instant…

As you may recall, I’m teaching analysis to this class of teachers, developing the $\epsilon$$\delta$ limit. Two weeks ago I bewildered everybody. Last week and this week, I set out to bewilder everyone even further.

Let me say what I’m going for here. The $\epsilon$$\delta$ limit is a notoriously difficult definition.1 How to scaffold my class to handle this difficulty? I am banking on the following strategy: make them need the definition. Make them unsatisfied with anything less. Continue poking holes in their current understanding, continue showing them inconsistencies between what they believe and the language they have to describe it, till they have no choice but to try to build something new. Then, let them try to build it. If they build the very thing I’m going for, rejoice. If they build something equally precise and powerful, rejoice. If they cannot build either (the most likely outcome, since the “right answer” took the world mathematical community 150 years to come up with), then it will still make powerful sense to them because it satisfactorily answers a question they were already engaged in trying to answer. That’s the plan anyway.

I will leave you with the two problem sets from the last class, and the readings and presentation from this one. I am very proud of the presentation. After that, I’ll write down one new thought for where to take this.

We engaged people’s attempts to define infinite decimals from the previous class, then abruptly shifted topics:

I let them work long enough so everyone got to do the first section of problems. My goals were:

1) Make participants recognize that they believe the speed of a moving object is something that exists in a particular moment of time.
2) Make them recognize that their naive definition of speed (distance / time) doesn’t actually handle this case.
3) Realize that we thus have a similar definitional problem as with repeating decimals.

We got this far. Then, with just 7 or so minutes left, I gave them another problem set:2

This problem set was designed to get somebody who has never studied calculus basically to take a simple derivative, to bring them into the conversation, and to refresh everyone else’s memory about the basic idea of derivatives. The last problem was on there just so that the calculus folks had a challenge available if they wanted it. Anyway, I had people finish the “Algebra Calisthenics” and “Speed” sections for homework.

This class, we began by engaging this homework, getting a feel for the standard calculus computation in which you identify the speed of an object in a moment as the value toward which average speeds seem to be headed as you look at smaller and smaller intervals. Then we began to press on what this really means.

I handed out a xerox of the scholium from the end of the first section of Book 1 of Newton’s Principia. (The last page of this pdf.) This is where Newton tries to explain what the hell he’s even talking about. I directed their attention to this telling sentence:

An in like manner, by the ultimate ratio of evanescent quantities is to be understood the ratio of the quantities, not before they vanish, nor afterwards, but with which they vanish.

Then, I showed them the following presentation. Wanting to share this with you is the real reason for this blog post. I had a lot of fun making it.

What’sCalculusReallyDoing (as pdf)

What’sCalculusReallyDoing (as powerpoint)

Then I passed out a choice excerpt from the awesome criticism of early calculus by Bishop George Berkeley. (Specif, section XIV.)

I asked for the connection between the definitional problem we have here and the definitional problem we had 2 classes ago regarding infinite decimals. (“They both involve getting closer and closer to something but never getting there.”) Then I asked them to try to come up with definitions to address these problems.

This is such a non-sequitur but here’s my one additional thought. I’ve been thinking about how to push participants to recognize a definition as unsatisfying. Tonight, reading Judith Grabiner’s 1983 essay in the AMM about Cauchy and the origins of the $\epsilon$$\delta$ limit (here it is as a pdf), I had an idea that is totally new to me. Retrospectively I think it’s sort of obvious, but I totally never thought of it before:

To get people to recognize that a definition is mathematically inadequate, have them try to use the definition, for example to prove something! In my case, all of them think that 1/3 = 0.333… Great. So, if we have a candidate definition of the meaning of limits or convergence, can we use it to prove 1/3 = 0.333…? If not, maybe we need a better definition.

(I had this idea when I read Grabiner’s statement that thought Cauchy gave the definition of the limit purely verbally and a bit vaguely, he translated it into the more rigorous language of inequalities when he actually started using it to prove theorems.)

[1] This is for at least 2 distinct (though related reasons): first of all, it’s got three nested quantifiers. “For all $\epsilon>0$, there exists a $\delta>0$, such that for all $x$ satisfying …” That just makes it inherently confusing. Secondly, it does not in any way psychologically resemble the intuitive image it is intended to capture. This is the definition of the limit. When I think of limits I have these beautiful visual images of little points getting closer to something. When I try to identify a limit, I just imagine the thing that they’re getting closer to. That’s the whole story. When I try to get rigorous, I replace this beautiful and simple image with three nested quantifiers. Yuck.

[2] You will notice some interconnections in the sequence of problems. After a few good experiences with this last year and then hearing how much fun everyone had at PCMI, I am beginning to feel like these sequences of densely but subtly interconnected problems are really, really awesome. Constructing them is a deep art and I am a tiny apprentice. But you can get started humbly and still see payoff: it was certainly a cool moment today in class when we went over these problems and a number of folks who had done out Speed problems #1-3 “the long way” realized that they could have applied their answer to Algebra Calisthenics #2 to do these three problems in moments in their heads.

I’m working with a new tutoring client, and therefore starting at the beginning in training her to stay engaged with the math rather than getting frustrated or trying to read the answer from my reactions to her guesses. The other night, I came up with a new metaphor to help her with this as she was trying to calculate the area of a circle section. I can’t believe I never thought of it before, it’s so obvious.

[N is visibly struggling to unify the geometric and algebraic information. I love it, I feel like I can literally see her brain growing, but she’s getting frustrated.]

Me: You’re going to grow from this.

N [skeptical]: Really?

Me: Yeah. You’re going to figure this out, and then you’re going to understand that you already had everything you need to figure it out. Have you seen Star Wars?

N: No, but I know it.

Me: You’re Luke, this is Vader. You face Vader and then you become a Jedi.

[Long pause while N thinks about the problem, punctuated by occasional exchanges like, “if I divide the circle area by $\theta$, does that give me the area of the wedge?” “I don’t know, make up numbers.”]

N: Oh! 360 divided by $\theta$ will give me the fraction of the circle that’s the wedge.

[4 or 5 second pause]

N: Right?

[N has evidently been watching my face intently for the last 4 or 5 seconds, trying to get external confirmation of her insight. When it isn’t forthcoming, she begins to doubt.]

Me: The reason I’ve been looking out the window is so you won’t try to get a cue off of my face. [Stands up] I’m going to do the more radical version of this and get out of your visual field. [Starts to head out of the room.]

N: You’re not going to tell me if I’m right?

Me: You’re trying to get me to fight Vader for you.

# 0.99999…

So I’m teaching this course this year. It’s for the math faculty of a high school. It’s called:

MA600 Algebra and Analysis with Connections to the K-12 Curriculum

I am unspeakably excited, and want to do the best job possible.

The class: 7 teachers, deeply committed to kids, serious, not real talkative, rightly protective about their time, which is in short supply, but eager to get sh*t done.

The content: Basically, all of mathematics, seen as a unified whole.

It’s met twice. The second class was last Thursday. I need to get my thoughts sorted out here. I’m expecting this to help me visualize the next moves more clearly, just by doing it, but I’d love your thoughts too.

I didn’t really know anything about the mathematical background of the group when I wrote the syllabus, so for the first class I gave them a getting-to-know-you problem set with a wide range of problems and just let them work the whole time. Magically the experience of watching folks work on the problems and then later looking at what they did on paper gave me just enough information to plan the direction of the class’ first unit. We’re beginning with analysis. My first goal: the $\epsilon$$\delta$ definition of the limit. (I.e., the definition of the limit, for the snobby among you.) My second: the completeness axiom.

The plan: generate the need to define the limit by working with 2 everyday concepts that are actually limits. Namely, infinite decimals, and instantaneous speed. My hope is that by pressing on these concepts, we’ll see that in spite of our familiarity with them, we don’t actually understand them unless we have a precise way to talk about limits. Then, develop the definition out of the need to fully understand the familiar. Then, develop the completeness axiom out of the desire to make sure infinite decimals have a limit.

Here’s what we did:

I opened class with a problem set designed to get them thinking about the meaning of decimals in particular, and various other contexts for the idea of limits. I shamelessly bit the format from PCMI. The problems span a wide range of skills and I didn’t leave enough time to do them all, so people could attack problems appropriate to their skill level. This is now my favorite way to differentiate problem sets, a propos of a) using it in some NYMC workshops last year, and b) hearing about how wonderful it was for everyone at PCMI.

Then, since we are all just getting to know each other, I did a short presentation on the mindset I wanted us to be in:

(Scribd did not handle slide 6 very well, which is too bad because I was proud of that slide. This is my first PowerPoint presentation ever. Actually I did it in Keynote.)

Then, we got to business. I put this up:

I asked them to talk about it with their tables. (I had them in 2 pairs and a group of three, in three tables in a horseshoe shape in front of the board. I like this and think I’ll keep it. Easy transitions from pair/group to whole-class; tables feel separate enough so you don’t feel like your conversation with your partner is in front of everybody; but everyone’s close enough so we can all talk. On day 1 I put us all around one table, for a sense of collegiality and common purpose, but it was too close; you couldn’t discreetly check in with your neighbor, for example.)

There was a widespread sense of mathematical discomfort, and rightly so. Infinite decimals enter most people’s math educations with no attention to the fact that they actually violate everything you’ve learned about math up to that point. You don’t get the full story until analysis, but unless you really get intimate with and own that content, you probably don’t connect what you learn there to what your teacher introduced without comment somewhere between 3rd and 7th grade, as though it weren’t a mind-boggling idea. “When you expand 1/3 as a decimal, the 3’s just keep going.” Or, “3.14159… It never ends or repeats.” Um, excuse me? It NEVER ENDS?

So it’s no surprise everybody has an underdeveloped idea of infinite decimals, and therefore that objects like 0.99999… cause some dissonance. This is very productive dissonance. I’m hoping it carries us all the way to the completeness axiom; we’ll see.

One of the three tables produced the standard argument that if x = 0.9999…, then 10x = 9.9999…, so 9x = 10x – x = 9, so x = 1; but even this table found this conclusion unsatisfying. I asked them why. The table that had produced the argument said, “usually this method gives you a fraction.” I asked for an example. They produced one from the problem set:

x = 1.363636…
100x = 136.363636…
99x = 100x – x = 135
x = 135/99 = 15/11

I asked how many folks found this argument convincing. 7 out of 7. (Well, one raised hand was kind of hesitant.)

I asked the same question about the same argument with .9999…. 4 out of 7. Then I dropped this:

How many people found this one convincing? 0 out of 7.

Reasonable.

Right?

Then what’s the difference?

At first, they cast about a bit, but then one of them said, “1.363636… has a finite limit, but …9999.0 doesn’t.” Their ideas began to coalesce around this type of language. Another one said, “we can actually estimate 1.363636…, for example we know it’s between 1 and 2.”

From the point of view I am ultimately heading for, this is the rub. Infinite decimals suggest convergent series, and the standard way to give them meaning as real numbers is that they are equal to the limit of the convergent series they suggest. …9999.0 suggests a wildly divergent series, so it cannot become a real number in the same way. (To bring home that convergence is the heart of the matter: there is an alternative way to define distance between numbers, the 10-adic metric, according to which it is actually …9999.0 that has the convergent series, and in this alternative system the above proof is valid and it actually does equal -1.) What I’d like us to do is a) define limits precisely; b) use this to prove that when a series has a limit, you can do the above type of manipulations to find it; c) try to prove that the series suggested by an infinite decimal always has a limit; d) realize that we can’t prove this without articulating the completeness axiom; e) articulate the axiom; and f) prove from the axiom that any infinite decimal has a real number limit. (Somewhere along the line, produce an $\epsilon$$\delta$ proof that 0.9999… = 1.) Now, how to orchestrate this…

For next time I told them to try to craft a definition of the meaning of an infinite decimal 0.abcd… I gave them a few minutes just before the end to discuss this with their groups. I’m expecting to learn a lot about their thinking from what they come up with, but I’m not counting on anyone to have a mathematically satisfying answer. I’ll be pleased if somebody does though.

As I think about next class, here’s what’s on my mind:

1) When we develop the $\epsilon$$\delta$ limit, what I’m going for is for this definition to feel like a satisfying relief. I know how easy it is for this definition instead to feel like a horrible monstrosity designed to oppress analysis students. I think what I have to do is keep them thinking about the reasons why anything less than this definition is too vague, which means I need to keep coming up with objects and problems that throw monkey wrenches into whatever more naive definitions they go for. (Of course, if they come up with something equally precise as the $\epsilon$$\delta$ limit but different, that would be amazing.) I feel like we’re off to a good start on this, but I want a fuller catalogue of head-scratchers (like …9999.0 = -1) to push the level of precision higher.

2) Relatedly, I sense a danger that the “real answers” will be unsatisfying because it’ll feel like “wait, I already said that.” For example, the participant who said that the difference between 1.3636… = 15/11 and …9999.0 = -1 is that “the first one has a finite limit”… I mean this is basically the answer. But it’s not based on a precise definition of limit yet, so it’s not what I want yet. I’m afraid of a “what was the big deal?” moment when we’ve got the real sh*t up there. I think the way to avoid this lies in that catalogue of head-scratchers I need to develop, so that nothing less than the real thing is satisfying. What do you think?

3) Where to go immediately next. Basically the question is: stick with decimals? Or change gears completely and press on the notion of instantaneous speed? Most (not all, I think) of the teachers have had a calculus course, but think at most 1 or 2 of them have internalized the philosophical lesson that instantaneous speed needs to be defined as a limit in order for us to even access it. I’m attracted to the idea of switching gears because I’m drawn to the connection between the disparate realms: two highly familiar, but totally different, objects – infinite decimals and speed in a moment – both getting pressed on to the point where you realize you never fully understood either one, and then you realize that the missing idea you need is the same thing in the two cases. (A precise way to talk about what number some varying quantity is “heading toward.”)

Actually as I write this out, it seems clear to me that switching gears is the way to go. I think it’ll give us a clearer understanding of what we’re missing with the decimals. Also, it’ll allow us to access all this rich historical stuff around the development of calculus. For example, maybe I’ll share with them some choice quotes from Bishop Berkeley’s The Analyst, to help articulate why the 18th century definition of the derivative was inadequate.

Anyway. Very excited about all this. Will definitely keep you posted.

# The Math Wizard

Okay, one more shout out.

My colleague Japheth Wood (News from the Math Wizard), with whom I’m delighted to be co-teaching a class for preservice teachers, is an awesome problem composer. He’s also been dipping into the sea of math and math ed blogging, one toe at a time. He’s finally got one whole foot in:

Check it.

In particular, check out that image.

Don’t move past it to the text (mine or Japheth’s) until you’ve sat with it long enough to absorb everything. If you find yourself with a mathematical question, don’t move on till you’ve tried to answer it.

Seriously; stop reading and go look.

There’s a lot of difference, pedagogically and content-wise, between this image and Dan’s boat-in-the-river video but there’s something very important and very exciting in common. Both manage to ask a very specific and mathematically rich question without whispering a single word. I think the natural current could be strengthened a bit by putting little venus flytraps on top of the square numbers, but that’s the only improvement I can think of.

So, questions for you:

a) What question, if any, does that image leave you with? (Am I right that there is a natural question you can’t help but have once you’ve absorbed the image?)
b) What are the features of the image that lead to the question? Given a mathematical question, how do we go about turning it into a wordless image that asks it?

FOLLOWUP (10/1/10):

Japheth and I passed out slips of paper to our class of preservice math teachers a week ago. On the slips were either Japheth’s original image, or tieandjeans’ modification. We asked them to write down a mathematical question that the paper provoked, and then try to answer it. We didn’t give them a ton of time. (Less than 10 min.) Interestingly, while I thought the sense of danger in the modification would make the gravitational pull toward our intended question (will the grasshopper manage to avoid all the squares?) greater, our students’ knowledge of Super Mario Bros was a distraction, because the up/down motion of the plants, and the question of Mario’s specific trajectory, became relevant considerations. (You can see below that one group, perhaps reading what we were going for, explicitly ruled out those considerations.) So I think the students that got the original grasshopper image actually gravitated toward the intended question more predictably. I still think the sense of danger would help, but maybe we just keep the grasshopper and add venus flytraps that appear static and aren’t close to the trajectory?

Anyway, here’s what they came up with. As you can see, for all of the above, the natural current is still pretty strong.

Grasshopper / Mario Problems

—=—=—=—=—=—=—=—=—=—=—=—

If the placement of fly traps continues, and Mario times his jumps so the trajectory never hits the plant, will he ever land on a fly trap and die?

$a^2 = 4n + 2$. No.

—=—=—=—=—=—=—=—=—=—=—=—

Does Jiminy land on a box?

$4k+2$, $k\in\{0,1,2,\ldots\}$
$k^2 = 4k+2$
$k^2 - 4k - 2=0$ —> $k = 2 \pm 2\sqrt{6}$

—=—=—=—=—=—=—=—=—=—=—=—

When does Mario’s jump not clear the venus fly trap?

—=—=—=—=—=—=—=—=—=—=—=—

Will the grasshopper land on one of the empty boxes (perfect squares)?

—=—=—=—=—=—=—=—=—=—=—=—

Hopping by 4. What is the mathematical formula to determine where Mario lands?

Mario will never land on a perfect square.

—=—=—=—=—=—=—=—=—=—=—=—

How many times will the grasshopper land between consecutive perfect squares?

For # between $(n-1)^2$, $n^2$, if odd then $(n-1)/2$, if even then $n/2$.

Will he ever land on a perfect square?

—=—=—=—=—=—=—=—=—=—=—=—

Can the grasshopper keep jumping without hitting the black box? If no, then when will he hit the black box?

—=—=—=—=—=—=—=—=—=—=—=—

Patterns
Grey boxes increase by +4
Black boxes increase by consecutive odd #s, +3, +5, +7, +9

Questions
Do you notice a pattern between the grey boxes? How about the black? Can you predict what # Mario will land on next?

The next plant will appear at 36. Mario will land at 26, then 30.

When will Mario land on a flower (perfect square)?

—=—=—=—=—=—=—=—=—=—=—=—

$2+4y = n^2$

Will Mario die and when?

—=—=—=—=—=—=—=—=—=—=—=—

When will Mario land on a plant?

Can he change the size of his jumps and still ensure he will not land on a plant?

—=—=—=—=—=—=—=—=—=—=—=—

Will the grasshopper ever land on a black square?

Can (4n-2) be a perfect square?

—=—=—=—=—=—=—=—=—=—=—=—

Will the grasshopper land on a boxed number (i.e. a perfect square)?

$4n - 2 = (y^2 + 0^2)$
Sum of squares must be a multiple of 4 or odd.
No.

(Ed. note: they’re misquoting a result they found the previous week. The result was about the difference of squares.)

—=—=—=—=—=—=—=—=—=—=—=—


text:
MathOverflow will be down for maintenance for approximately 3 hours, starting Monday evening (06/24/2013) at approximately 9:00 PM Eastern time (UTC-4).

3 added 743 characters in body

Bochner's theorem states that a positive definite function is the Fourier transform of a finite Borel measure. As well, an easy converse of this is that a Fourier transform must be positive definite.

My question is: is there a high-brow explanation for why positive definiteness and Fourier transforms go hand-in-hand?

As I understand it, positive definiteness imposes wonderfully strong regularity conditions on the function. We immediately deduce that the function is bounded above at its value at 0, that it is non-negative at 0 and that continuity at 0 implies continuity everywhere.

A leading example I have in mind comes from probability. One can show (Levy's Theorem) that a sum of iid rv converges weakly to some probability distribution by considering the product of characteristic functions and showing that its tail converges to 1 around an interval containing 0, so by positive definiteness and by the identity $1-\mbox{Re} \phi(2t) \leq 4(1-\mbox{Re} \phi(t))$ this implies convergence to a degenerate distribution. It just seems rather mysterious to me how this kind of local regularity becomes global.

Edit:

To be a little more specific, I understand that the Radon Nikodym derivative is positive and $e^{ix}$ is positive definite. I am more interested in consequences of positive-definiteness on the regularity of the function. For example, if one takes the 2x2 positive definite matrix associated with the function and considers its determinant, it follows that $|f(x)|\leq |f(0)|$. If I take the 3x3 positive definite matrix, I can conclude that if $f$ is continuous at 0, it is then continuous everywhere. My issue is that these types of arguments give me no intuition at all as to what positive definiteness is.

2 edited title

1

# Demystifying Positive Definite Functions

Bochner's theorem states that a positive definite function is the Fourier transform of a finite Borel measure. As well, an easy converse of this is that a Fourier transform must be positive definite.

My question is: is there a high-brow explanation for why positive definiteness and Fourier transforms go hand-in-hand?

As I understand it, positive definiteness imposes wonderfully strong regularity conditions on the function. We immediately deduce that the function is bounded above at its value at 0, that it is non-negative at 0 and that continuity at 0 implies continuity everywhere.

A leading example I have in mind comes from probability. One can show (Levy's Theorem) that a sum of iid rv converges weakly to some probability distribution by considering the product of characteristic functions and showing that its tail converges to 1 around an interval containing 0, so by positive definiteness and by the identity $1-\mbox{Re} \phi(2t) \leq 4(1-\mbox{Re} \phi(t))$ this implies convergence to a degenerate distribution. It just seems rather mysterious to me how this kind of local regularity becomes global.