|
5 | 5 | AUTHORS:
|
6 | 6 |
|
7 | 7 | - Kevin Stueve (2010-01-17): in ``is_prime(n)``, delegated calculation to ``n.is_prime()``
|
8 |
| -- Lorenz Panny (2022): :class:`ProductTree`, :func:`prod_with_derivative` |
9 | 8 | """
|
10 | 9 |
|
11 | 10 | # ****************************************************************************
|
@@ -6286,249 +6285,3 @@ def dedekind_psi(N):
|
6286 | 6285 | """
|
6287 | 6286 | N = Integer(N)
|
6288 | 6287 | return Integer(N * prod(1 + 1 / p for p in N.prime_divisors()))
|
6289 |
| - |
6290 |
| - |
6291 |
| -class ProductTree: |
6292 |
| - r""" |
6293 |
| - A simple binary product tree, i.e., a tree of ring elements in |
6294 |
| - which every node equals the product of its children. |
6295 |
| - (In particular, the *root* equals the product of all *leaves*.) |
6296 |
| -
|
6297 |
| - Product trees are a very useful building block for fast computer |
6298 |
| - algebra. For example, a quasilinear-time Discrete Fourier Transform |
6299 |
| - (the famous *Fast* Fourier Transform) can be implemented as follows |
6300 |
| - using the :meth:`remainders` method of this class:: |
6301 |
| -
|
6302 |
| - sage: F = GF(65537) |
6303 |
| - sage: a = F(1111) |
6304 |
| - sage: assert a.multiplicative_order() == 1024 |
6305 |
| - sage: R.<x> = F[] |
6306 |
| - sage: ms = [x - a^i for i in range(1024)] # roots of unity |
6307 |
| - sage: ys = [F.random_element() for _ in range(1024)] # input vector |
6308 |
| - sage: zs = ProductTree(ms).remainders(R(ys)) # compute FFT! |
6309 |
| -
|
6310 |
| - This class encodes the tree as *layers*: Layer `0` is just a tuple |
6311 |
| - of the leaves. Layer `i+1` is obtained from layer `i` by replacing |
6312 |
| - each pair of two adjacent elements by their product, starting from |
6313 |
| - the left. (If the length is odd, the unpaired element at the end is |
6314 |
| - simply copied as is.) This iteration stops as soon as it yields a |
6315 |
| - layer containing only a single element (the root). |
6316 |
| -
|
6317 |
| - .. NOTE:: |
6318 |
| -
|
6319 |
| - Use this class if you need the :meth:`remainders` method. |
6320 |
| - To compute just the product, :func:`prod` is likely faster. |
6321 |
| -
|
6322 |
| - INPUT: |
6323 |
| -
|
6324 |
| - - ``leaves`` -- a sequence of elements in a common ring |
6325 |
| -
|
6326 |
| - EXAMPLES:: |
6327 |
| -
|
6328 |
| - sage: from sage.arith.misc import ProductTree |
6329 |
| - sage: R.<x> = GF(101)[] |
6330 |
| - sage: vs = [x - i for i in range(1,10)] |
6331 |
| - sage: tree = ProductTree(vs) |
6332 |
| - sage: tree.root() |
6333 |
| - x^9 + 56*x^8 + 62*x^7 + 44*x^6 + 47*x^5 + 42*x^4 + 15*x^3 + 11*x^2 + 12*x + 13 |
6334 |
| - sage: tree.remainders(x^7 + x + 1) |
6335 |
| - [3, 30, 70, 27, 58, 72, 98, 98, 23] |
6336 |
| - sage: tree.remainders(x^100) |
6337 |
| - [1, 1, 1, 1, 1, 1, 1, 1, 1] |
6338 |
| -
|
6339 |
| - :: |
6340 |
| -
|
6341 |
| - sage: vs = prime_range(100) |
6342 |
| - sage: tree = ProductTree(vs) |
6343 |
| - sage: tree.root().factor() |
6344 |
| - 2 * 3 * 5 * 7 * 11 * 13 * 17 * 19 * 23 * 29 * 31 * 37 * 41 * 43 * 47 * 53 * 59 * 61 * 67 * 71 * 73 * 79 * 83 * 89 * 97 |
6345 |
| - sage: tree.remainders(3599) |
6346 |
| - [1, 2, 4, 1, 2, 11, 12, 8, 11, 3, 3, 10, 32, 30, 27, 48, 0, 0, 48, 49, 22, 44, 30, 39, 10] |
6347 |
| -
|
6348 |
| - We can access the individual layers of the tree:: |
6349 |
| -
|
6350 |
| - sage: tree.layers |
6351 |
| - [(2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97), |
6352 |
| - (6, 35, 143, 323, 667, 1147, 1763, 2491, 3599, 4757, 5767, 7387, 97), |
6353 |
| - (210, 46189, 765049, 4391633, 17120443, 42600829, 97), |
6354 |
| - (9699690, 3359814435017, 729345064647247, 97), |
6355 |
| - (32589158477190044730, 70746471270782959), |
6356 |
| - (2305567963945518424753102147331756070,)] |
6357 |
| - """ |
6358 |
| - def __init__(self, leaves): |
6359 |
| - r""" |
6360 |
| - Initialize a product tree having the given ring elements |
6361 |
| - as its leaves. |
6362 |
| -
|
6363 |
| - EXAMPLES:: |
6364 |
| -
|
6365 |
| - sage: from sage.arith.misc import ProductTree |
6366 |
| - sage: vs = prime_range(100) |
6367 |
| - sage: tree = ProductTree(vs) |
6368 |
| - """ |
6369 |
| - V = tuple(leaves) |
6370 |
| - self.layers = [V] |
6371 |
| - while len(V) > 1: |
6372 |
| - V = tuple(prod(V[i:i+2]) for i in range(0,len(V),2)) |
6373 |
| - self.layers.append(V) |
6374 |
| - |
6375 |
| - def __len__(self): |
6376 |
| - r""" |
6377 |
| - Return the number of leaves of this product tree. |
6378 |
| -
|
6379 |
| - EXAMPLES:: |
6380 |
| -
|
6381 |
| - sage: from sage.arith.misc import ProductTree |
6382 |
| - sage: R.<x> = GF(101)[] |
6383 |
| - sage: vs = [x - i for i in range(1,10)] |
6384 |
| - sage: tree = ProductTree(vs) |
6385 |
| - sage: len(tree) |
6386 |
| - 9 |
6387 |
| - sage: len(tree) == len(vs) |
6388 |
| - True |
6389 |
| - sage: len(tree.remainders(x^2)) |
6390 |
| - 9 |
6391 |
| - """ |
6392 |
| - return len(self.layers[0]) |
6393 |
| - |
6394 |
| - def __iter__(self): |
6395 |
| - r""" |
6396 |
| - Return an iterator over the leaves of this product tree. |
6397 |
| -
|
6398 |
| - EXAMPLES:: |
6399 |
| -
|
6400 |
| - sage: from sage.arith.misc import ProductTree |
6401 |
| - sage: R.<x> = GF(101)[] |
6402 |
| - sage: vs = [x - i for i in range(1,10)] |
6403 |
| - sage: tree = ProductTree(vs) |
6404 |
| - sage: next(iter(tree)) == vs[0] |
6405 |
| - True |
6406 |
| - sage: list(tree) == vs |
6407 |
| - True |
6408 |
| - """ |
6409 |
| - return iter(self.layers[0]) |
6410 |
| - |
6411 |
| - def root(self): |
6412 |
| - r""" |
6413 |
| - Return the value at the root of this product tree (i.e., the product of all leaves). |
6414 |
| -
|
6415 |
| - EXAMPLES:: |
6416 |
| -
|
6417 |
| - sage: from sage.arith.misc import ProductTree |
6418 |
| - sage: R.<x> = GF(101)[] |
6419 |
| - sage: vs = [x - i for i in range(1,10)] |
6420 |
| - sage: tree = ProductTree(vs) |
6421 |
| - sage: tree.root() |
6422 |
| - x^9 + 56*x^8 + 62*x^7 + 44*x^6 + 47*x^5 + 42*x^4 + 15*x^3 + 11*x^2 + 12*x + 13 |
6423 |
| - sage: tree.root() == prod(vs) |
6424 |
| - True |
6425 |
| - """ |
6426 |
| - assert len(self.layers[-1]) == 1 |
6427 |
| - return self.layers[-1][0] |
6428 |
| - |
6429 |
| - def leaves(self): |
6430 |
| - r""" |
6431 |
| - Return a tuple containing the leaves of this product tree. |
6432 |
| -
|
6433 |
| - EXAMPLES:: |
6434 |
| -
|
6435 |
| - sage: from sage.arith.misc import ProductTree |
6436 |
| - sage: R.<x> = GF(101)[] |
6437 |
| - sage: vs = [x - i for i in range(1,10)] |
6438 |
| - sage: tree = ProductTree(vs) |
6439 |
| - sage: tree.leaves() |
6440 |
| - (x + 100, x + 99, x + 98, ..., x + 93, x + 92) |
6441 |
| - sage: tree.leaves() == tuple(vs) |
6442 |
| - True |
6443 |
| - """ |
6444 |
| - return self.layers[0] |
6445 |
| - |
6446 |
| - def remainders(self, x): |
6447 |
| - r""" |
6448 |
| - Given a value `x`, return a list of all remainders of `x` |
6449 |
| - modulo the leaves of this product tree. |
6450 |
| -
|
6451 |
| - The base ring must support the ``%`` operator for this |
6452 |
| - method to work. |
6453 |
| -
|
6454 |
| - INPUT: |
6455 |
| -
|
6456 |
| - - ``x`` -- an element of the base ring of this product tree |
6457 |
| -
|
6458 |
| - EXAMPLES:: |
6459 |
| -
|
6460 |
| - sage: from sage.arith.misc import ProductTree |
6461 |
| - sage: vs = prime_range(100) |
6462 |
| - sage: tree = ProductTree(vs) |
6463 |
| - sage: n = 1085749272377676749812331719267 |
6464 |
| - sage: tree.remainders(n) |
6465 |
| - [1, 1, 2, 1, 9, 1, 7, 15, 8, 20, 15, 6, 27, 11, 2, 6, 0, 25, 49, 5, 51, 4, 19, 74, 13] |
6466 |
| - sage: [n % v for v in vs] |
6467 |
| - [1, 1, 2, 1, 9, 1, 7, 15, 8, 20, 15, 6, 27, 11, 2, 6, 0, 25, 49, 5, 51, 4, 19, 74, 13] |
6468 |
| - """ |
6469 |
| - X = [x] |
6470 |
| - for V in reversed(self.layers): |
6471 |
| - X = [X[i//2] % V[i] for i in range(len(V))] |
6472 |
| - return X |
6473 |
| - |
6474 |
| - |
6475 |
| -def prod_with_derivative(pairs): |
6476 |
| - r""" |
6477 |
| - Given a list of pairs `(f, \partial f)` of ring elements, return |
6478 |
| - the pair `(\prod f, \partial \prod f)`, assuming `\partial` is an |
6479 |
| - operator obeying the standard product rule. |
6480 |
| -
|
6481 |
| - This function is entirely algebraic, hence still works when the |
6482 |
| - elements `f` and `\partial f` are all passed through some ring |
6483 |
| - homomorphism first. One particularly useful instance of this is |
6484 |
| - evaluating the derivative of a product of polynomials at a point |
6485 |
| - without fully expanding the product; see the second example below. |
6486 |
| -
|
6487 |
| - INPUT: |
6488 |
| -
|
6489 |
| - - ``pairs`` -- a sequence of tuples `(f, \partial f)` of elements |
6490 |
| - of a common ring |
6491 |
| -
|
6492 |
| - ALGORITHM: Repeated application of the product rule. |
6493 |
| -
|
6494 |
| - EXAMPLES:: |
6495 |
| -
|
6496 |
| - sage: from sage.arith.misc import prod_with_derivative |
6497 |
| - sage: R.<x> = ZZ[] |
6498 |
| - sage: fs = [x^2 + 2*x + 3, 4*x + 5, 6*x^7 + 8*x + 9] |
6499 |
| - sage: prod(fs) |
6500 |
| - 24*x^10 + 78*x^9 + 132*x^8 + 90*x^7 + 32*x^4 + 140*x^3 + 293*x^2 + 318*x + 135 |
6501 |
| - sage: prod(fs).derivative() |
6502 |
| - 240*x^9 + 702*x^8 + 1056*x^7 + 630*x^6 + 128*x^3 + 420*x^2 + 586*x + 318 |
6503 |
| - sage: F, dF = prod_with_derivative((f, f.derivative()) for f in fs) |
6504 |
| - sage: F |
6505 |
| - 24*x^10 + 78*x^9 + 132*x^8 + 90*x^7 + 32*x^4 + 140*x^3 + 293*x^2 + 318*x + 135 |
6506 |
| - sage: dF |
6507 |
| - 240*x^9 + 702*x^8 + 1056*x^7 + 630*x^6 + 128*x^3 + 420*x^2 + 586*x + 318 |
6508 |
| -
|
6509 |
| - The main reason for this function to exist is that it allows us to |
6510 |
| - *evaluate* the derivative of a product of polynomials at a point |
6511 |
| - `\alpha` without ever fully expanding the product *as a polynomial*:: |
6512 |
| -
|
6513 |
| - sage: alpha = 42 |
6514 |
| - sage: F(alpha) |
6515 |
| - 442943981574522759 |
6516 |
| - sage: dF(alpha) |
6517 |
| - 104645261461514994 |
6518 |
| - sage: us = [f(alpha) for f in fs] |
6519 |
| - sage: vs = [f.derivative()(alpha) for f in fs] |
6520 |
| - sage: prod_with_derivative(zip(us, vs)) |
6521 |
| - (442943981574522759, 104645261461514994) |
6522 |
| - """ |
6523 |
| - class _aux: |
6524 |
| - def __init__(self, f, df): |
6525 |
| - self.f, self.df = f, df |
6526 |
| - |
6527 |
| - def __mul__(self, other): |
6528 |
| - return _aux(self.f * other.f, self.df * other.f + self.f * other.df) |
6529 |
| - |
6530 |
| - def __iter__(self): |
6531 |
| - yield self.f |
6532 |
| - yield self.df |
6533 |
| - |
6534 |
| - return tuple(prod(_aux(*tup) for tup in pairs)) |
0 commit comments