Skip to content

Instantly share code, notes, and snippets.

@pjt33
pjt33 / linear_diophantine.py
Created Jan 21, 2020
Count solutions to linear Diophantine equations
View linear_diophantine.py
from collections import Counter
from fractions import Fraction
def _gcd(a, b):
while a:
a, b = b % a, a
return b
@pjt33
pjt33 / math3444274.py
Last active Nov 23, 2019
Dragon merging
View math3444274.py
#!/usr/bin/python3
# https://math.stackexchange.com/q/3444274
from fractions import Fraction
def build_states(target_dragon):
states = {}
def toBinary(x, width):
View math3414931.py
#!/usr/bin/python3
from functools import lru_cache
@lru_cache(None)
def binom(n, k):
if n < 0 or k < 0 or k > n:
return 0
rv = 1
for i in range(k):
@pjt33
pjt33 / CR228847.numpy
Last active Sep 13, 2019
Chebyshev pseudoprime test
View CR228847.numpy
#!/usr/bin/python3
import timeit
from numpy import mod
from numpy.polynomial.polynomial import Polynomial
def is_prime_naive(n):
if n <= 2 or n % 2 == 0:
return n == 2
@pjt33
pjt33 / math3343318.md
Last active Sep 3, 2019
Strings with prefix imbalance
View math3343318.md

Let's define $a_{n,m}$ to be the number of words of length $n$ with a surplus of $m$ A's. E.g. with $k=2$, AAAB would contribute $1$ to $a_{4,1}$; with $k=3$ it would contribute $1$ to $a_{4,0}$.

Given a word of length $n-1$ we can always append an A to get a valid word; we can append a B to get a valid word iff the surplus is at least $k$. So if we set up the generating function $$f(x,z) = \sum_{n \ge 0} \sum_{m=0}^n a_{n,m} x^n z^m$$ we find that $$f(x, z) = 1 + xz f(x, z) + xz^{-p} \sum_{n \ge 0} \sum_{m=p}^n a_{n,m} x^n z^m \tag{1}$$ which doesn't look very promising.

However, I believe Deutsch proved that this kind of regression always gives a Riordan array, so let's assume that we can write $f(x, z)$ in the form $$f(x, z) = \sum_{m \ge 0} z^m d(x) (xh(x))^m \tag{2}$$

Useful observation: the length minus the surplus is always a multiple of $k+1$ (easily shown by induction), so $d(x)$ only has non-zero coefficients at powers of $x^{k+1}$, or $d(x) = d_k(x^{k+1})$.

Useful observation: if we ext

@pjt33
pjt33 / DriveYaNuts.cs
Created Sep 3, 2019
Drive Ya Nuts solution counter
View DriveYaNuts.cs
using System;
using System.Collections.Generic;
using System.Linq;
// WARNING: THIS IS A QUICK HACK, NOT PRODUCTION CODE. IT CAN PROBABLY BE OPTIMISED QUITE SIGNIFICANTLY.
namespace Sandbox
{
// https://math.stackexchange.com/q/3306917
static class DriveYaNuts
@pjt33
pjt33 / MO338802.cs
Created Aug 26, 2019
Perfectly nonlinear involutions
View MO338802.cs
using System;
using System.Collections.Generic;
using System.Linq;
using System.Numerics;
namespace Sandbox
{
class MO338802
{
public static void Main()
@pjt33
pjt33 / Carpsen.swift
Created Jan 14, 2019
Eratospheres in Swift
View Carpsen.swift
func eratosthenesSieve(to n: Int) -> [Int] {
guard 2 <= n else { return [] }
var composites = Array(repeating: false, count: n + 1)
var primes: [Int] = []
let d = Double(n)
let upperBound = Int((d / _log(d)) * (1.0 + 1.2762/_log(d)))
primes.reserveCapacity(upperBound)
let squareRootN = Int(d.squareRoot())
View mathse3015816.py
#!/usr/bin/python3
import math
def phi(n):
# Naïve prime factorisation
ps = set()
m = n
while m % 2 == 0:
@pjt33
pjt33 / math2962644.sage
Last active Oct 20, 2018
Number of real zeroes of iterated polynomial: x^3-2x+1
View math2962644.sage
from sage.rings.polynomial.real_roots import *
x = polygen(ZZ)
def conv_p(p, n):
if n == 0:
return x
else:
return p(conv_p(p, n-1))
for n in range(1, 21): print add(map(lambda (a,b): b, real_roots(conv_p(x^3-2*x-1, n))))
You can’t perform that action at this time.