Godtone

Joined 17 December 2020

User page
Discussion

← Older edit Newer edit →

VisualWikitext

Retrieved from "https://en.xen.wiki/w/User:Godtone"

@@ Line 381: / Line 381: @@
 |}
 <references/>
+== My Python 3 code ==
+<syntaxhighlight lang="python">
+import math
+import functools
+# prod() computes the product of a list of things.
+# It can be thought of as the multiplicative counterpart to sum().
+# EG: `prod([3,2,4])` evaluates to (3*2)*4 = 24.
+prod = functools.partial(functools.reduce, lambda x,y: x*y)
+# p is a complete list of (at least 2) primes up to and including some prime.
+# gen_prime(p) appends the next prime to p and returns that prime.
+def gen_prime(p):
+	c = p[-1] + 2
+	while True:
+		i = 0
+		found_divisor = False
+		while c >= p[i]*p[i]:
+			if c % p[i] == 0:
+				found_divisor = True
+				break
+			i+=1
+		if not found_divisor: # aka if c is prime
+			p.append(c)
+			return c
+		c += 2
+primes = [2,3]
+# Cannot give index error (unless you run out of memory) and
+# allows implicit convenient access to global variable `primes`
+def prime_idx(i):
+	global primes
+	while i >= len(primes):
+		gen_prime(primes)
+	return primes[i]
+# Returns a list of primes up to and (if applicable) including the limit,
+# representing the corresponding prime limit.
+def prime_limit(limit): # includes limit (if applicable)
+	l = []
+	i = 0
+	while prime_idx(i) <= limit:
+		l.append(prime_idx(i))
+		i += 1
+	return l
+# Can be thought of as the conceptual opposite of strip_list_of_right().
+def right_justify_list(l, n, justify_with=0):
+	while len(l) < n:
+		l.append(justify_with)
+	return l
+# Can be thought of as the conceptual opposite of right_justify_list().
+# Note that what_to_strip is a list of elements eligible for stripping (removal).
+def strip_list_of_right(l, what_to_strip=[0]):
+	while len(l)>0 and (l[-1] in what_to_strip):
+		l.pop()
+	return l
+# If p is None (default):
+#    Factorises a positive integer into a list of the exponents/powers
+#    of corresponding primes starting with the exponent/power of 2.
+#    Note that any non-positive integer will give an empty list/factorisation
+#    and that this is the same factorisation that the number 1 gets.
+# Otherwise:
+#    Attempts factorisation based on a list of positive integers p
+#    which can be used to try to divide n down to 1.
+#    If n cannot be divided down to 1, an empty list is returned.
+#    Note that two versions of each function dealing with factorisations exist;
+#    one of which specifying such a list of positive integers p.
+#    Also note that the positive integers are tried in the order they are listed.
+def fact_int(n, p=None):
+	if p==None:
+		f = []
+		while n > 1:
+			f.append(0)
+			while n % prime_idx(len(f)-1) == 0:
+				f[-1] += 1
+				n //= prime_idx(len(f)-1)
+		return f
+	else:
+		f = [0]*len(p)
+		i = 0
+		while n > 1 and i < len(p):
+			while n % p[i] == 0:
+				n //= p[i]
+				f[i] += 1
+			i += 1
+		return f if n==1 else []
+# Note that this works for negative exponents of primes too but produces
+# a floating point value rather than an exact rational representation.
+def unfact_int(f, p=None):
+	return (
+if f==[]
+		else prod([ prime_idx(i)**f[i] for i in range(len(f)) ]) if p==None
+		else prod([ p[i]**f[i] for i in range(len(f)) ])
+		)
+# Note: "fact" stands for "factorise(d form (of))/factorisation".
+# Note: when combining two factorisations in some way,
+#       both must either be normal or subgroup-based.
+def mul_fact(x, y):
+	x, y = x.copy(), y.copy()
+	max_len = max(len(x),len(y))
+	x, y = right_justify_list(x,max_len), right_justify_list(y,max_len)
+	result = [x[i]+y[i] for i in range(max_len)]
+	return result
+def recip_fact(f, p=None): # useless parameter p used for niceness
+	return [-i for i in f]
+def div_fact(x, y):
+	return mul_fact(x, recip_fact(y))
+def fact(r, p=None):
+	return (
+		strip_list_of_right(div_fact( fact_int(r[0]), fact_int(r[1]) ))
+		if p==None else
+		div_fact( fact_int(r[0],p), fact_int(r[1],p) )
+		)
+def unfact(f, p=None):
+	return unfact_int([max(n,0) for n in f],p), unfact_int([-min(d,0) for d in f],p)
+def iv(n, d, p=None): # provides automatic reduction to simplest form
+	return unfact(fact( (n,d),p ),p)
+def as_float(r,p=None):
+	if type(r)==list:
+		uf = unfact(r, p)
+		return uf[0] / uf[1]
+	elif type(r)==tuple:
+		return r[0] / r[1]
+	else:
+		return r
+# Converts a prime factorisation into a prettier and potentially
+# more readable (simple) string representation. (Not for subgroup-factorisations.)
+# redundancy is either 0, 1 or 2 and has the following effects:
+# redundancy=0 (default):
+#    Only includes primes with nonzero exponent and omits the exponent if
+#    if the exponent is 1 (AKA if the prime would be followed by '^1').
+# redundancy=1:
+#    Same as redundancy=0 except exponents are always explicit ('^1' is included).
+# redundancy=2:
+#    Shows all primes up to the largest prime with nonzero exponent.
+def fact_to_str(f, redundancy=0):
+	sl = []
+	for i in range(len(f)):
+		if redundancy==2 or f[i]!=0:
+			sl.append( str(prime_idx(i))
+				+( '^'+str(f[i]) if redundancy>0 or f[i]!=1 else '' )
+				)
+	if len(sl) == 0:
+		return '1'
+	return ' * '.join(sl)
+# The closest approximation of r in (d)ED(n) is `round(in_ed(r,d,n))`, or
+# more commonly, the closest approximation of r in n-EDO is `round(in_ed(r,n))`.
+def steps(r, et2):
+	return math.log(r,2) / et2
+# Note that error is measured as deviation from the exact value of r in
+# (divisions)ED(octave_equivalent). The unsigned AKA absolute error of the
+# closest approximation of r in n-EDO would be written `abs(err_in_ed(r,n))`.
+def step_err(r, et2):
+	exact = steps(r, et2)
+	return round(exact) - exact
+# For use with steps(), step_err(), etc. as the et2 argument
+def ed(equal_divisions, of_n=2):
+	return math.log(of_n,2)/equal_divisions
+def pval(r, et2, p=None):
+	if p==None:
+		p = prime_limit(prime_idx( len(r)-1 ))
+	return sum([ r[i]*round(steps(p[i],et2)) for i in range(len(r)) ])
+</syntaxhighlight>