Hace unos días [ver nota], Terence Tao propuso un reto para descomponer 300K! como el producto de 300K factores mayores que 100K. Un ejemplo más pequeño es descomponer 10! como el producto de 10 factores mayores o iguales que 3.
10! = 1 * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10
10! = 3 * 3 * 4 * 4 * 4 * 5 * 5 * 6 * 6 * 7
Él logró demostrar que hay una descomposición de 300K! como el producto de 300K factores mayores que solo 90K y propuso un método para intentar usar números más grandes. Esto es parte de un intento para demostrar una conjetura, por lo que el problema completo es más general, pero intentemos resolver este caso.
[Edición: Después de publicar esto en HN, en un comentario de btilly me mostró que Andrew Sutherland resolvió este caso un día después de su publicación. Pero usando un método muy diferente.]
[Notar: La versión en inglés fue publicada en abril, la versión en castellano en diciembre, pero la backdate para que quede bien en el calendario de posts.]
Fijaremos N=300K, que es el caso analizado por Tao en su post. La idea de Tao es comenzar con un número impar B, que es producto de 300K números impares mayores que 100K, e intentar corregir la diferencia. Un punto importante es que 300! es par, así que la idea es usar los 2s para ayudarnos en esta tarea.
Él define los primos B-heavy, que son los que aparecen más veces en B que en N!=300!, y los primos N!-heavy, que son los que aparecen más veces en N!=300! que en B. La idea es reemplazar cada primo B-heavy b por un primo N!-heavy n multiplicado por una potencia de 2, o directamente por una potencia pura de 2. Esto aumentará el factor de B y obtendremos la descomposición que queremos. Es importante usar la menor cantidad posible de 2 para no quedarnos sin ninguno. Todo esto está explicado mejor en el post de Tao, así que vale la pena leerla.
Él pudo aplicar este método para obtener una descomposición con números mayores de 90K, pero el apareamiento y el uso de 2 no fueron óptimos, por lo que no fue suficiente para 100K. Vamos a encontrar la manera óptima de realizar este apareamiento y obtener una descomposición con números mayores de 100K.
Preámbulo
Usaremos Racket. Primero, intentaremos reproducir el resultado con 90K. Fijaremos N=300K y escribiremos una versión memorizada de factorizar, porque la factorización es difícil y lenta, y reutilizaremos mucho los mismo números.
#lang racket
(require math/number-theory)
(define N 300000)
(define secret-factorize-cache (make-hash))
(define (factorize/cache n)
(hash-ref! secret-factorize-cache n (lambda () (factorize n))))
Factorización de N!=300K!
Es mucho mejor factorizar cada número y juntar las factorizaciones que usar el factorial interno de Racket e intentar factorizar el resultado, lo que podría llevar mucho tiempo (o fallar, porque factorizar números grandes es difícil).
(define (hash-add-factorization! h f)
(for ([pair (in-list f)])
(define p (car pair))
(define m (cadr pair))
(hash-update! h p {lambda (x) (+ x m)} 0)))
(define (hash->factorization h)
(define l (for/list ([(p m) (in-hash h)])
(list p m)))
(sort l < #:key car))
(define (factorize-factorial N)
(define h (make-hash))
(for ([n (in-range N)])
(hash-add-factorization! h (factorize/cache (add1 n))))
(hash->factorization h))
(define (factorization-total N)
(define t 0)
(for ([pm (in-list N)])
(set! t (+ t (cadr pm))))
t)
(display "N! factorization calculation (N=300000): ")
(define Nf2-factorization (time (factorize-factorial N)))
(display "Number of primes including 2: ")
(displayln (length Nf2-factorization))
(display "Number of primes with multiplicity including 2: ")
(displayln (factorization-total Nf2-factorization))
(display "Multiplicity of 2: ")
(define Nf-2multiplicity (if (and (not (null? Nf2-factorization))
(= (caar Nf2-factorization) 2))
(cadar Nf2-factorization)
0))
(displayln Nf-2multiplicity)
(display "Should be 299992. Difference: ")
(displayln (- 299992 Nf-2multiplicity))
(displayln "N! factorization excluding 2: ")
(define Nf-factorization (if (and (not (null? Nf2-factorization))
(= (caar Nf2-factorization) 2))
(cdr Nf2-factorization)
Nf2-factorization))
(display "Number of primes excluding 2: ")
(displayln (length Nf-factorization))
(display "Number of primes with multiplicity excluding 2: ")
(displayln (factorization-total Nf-factorization))
(display "Partial list excluding 2: ")
(displayln (append (take Nf-factorization 7)
(list "...")
(take-right Nf-factorization 3)))
(newline)
La salida es
N! factorization calculation (N=300000): cpu time: 4828 real time: 4825 gc time: 406
Number of primes including 2: 25997
Number of primes with multiplicity including 2: 1059501
Multiplicity of 2: 299992
Should be 299992. Difference: 0
N! factorization excluding 2:
Number of primes excluding 2: 25996
Number of primes with multiplicity excluding 2: 759509
Partial list excluding 2: ((3 149995) (5 74998) (7 49996) (11 29997) (13 24997) (17 18749) (19 16665) ... (299977 1) (299983 1) (299993 1))
Calcular la factorización tarda unos 4 segundos, y obtenemos el mismo número de 2s que Tao. (Muchos 2s).
Factorización de B (L=90K, A=50)
Para reproducir el resultado, veremos los mismos L y A que Tao, por lo que el número B son 50 copias de los números impares entre 90K y 102K. Mi cálculo considera también el caso en el que el número A no divide a N, y decidí agregar más números impares a partir de L.
(define (factorize-B-number N L size)
(define h (make-hash))
(define from (if (odd? L) L (add1 L)))
(for ([n (in-range N)])
(hash-add-factorization! h (factorize/cache (+ from (* (remainder n size) 2)))))
(hash->factorization h))
(display "B factorization calculation (L=90000, A=50): ")
(define B-factorization (time (factorize-B-number N 90000 (quotient N 50))))
(display "Number of primes: ")
(displayln (length B-factorization))
(display "Number of primes with multiplicity: ")
(displayln (factorization-total B-factorization))
(display "Partial list: ")
(displayln (append (take B-factorization 7)
(list "...")
(take-right B-factorization 3)))
(newline)
(display "B has more primes than N!, with multiplicity excluding 2: ")
(displayln (- (factorization-total B-factorization)
(factorization-total Nf-factorization)))
(display "Should be 14891. Difference: ")
(displayln (- (- (factorization-total B-factorization)
(factorization-total Nf-factorization))
14891))
(newline)
La salida es
B factorization calculation (L=90000, A=50): cpu time: 531 real time: 534 gc time: 46
Number of primes: 3099
Number of primes with multiplicity: 774450
Partial list: ((3 150000) (5 75000) (7 50000) (11 29900) (13 25000) (17 18700) (19 16600) ... (101977 50) (101987 50) (101999 50))
B has more primes than N!, with multiplicity excluding 2: 14941
Should be 14891. Difference: 50
Es más rápido y solo toma medio segundo porque las factorizaciones están memorizadas. Obtuve 50 primos de más, contados con multiplicidad. Estaba usando un método complicado para calcular los números impares que se repiten, pero lo cambié al método simple y lento. Sin embargo, los 50 primos adicionales no desaparecen. Agradecería cualquier sugerencia o reporte de errores.
N!-heavy y B-heavy
Ahora comparamos la factorización de N! y B, eliminando los primos repetidos (contados con multiplicidad).
(define (factorization-split N B)
(define h (make-hash))
(for ([pm (in-list N)])
(hash-update! h (car pm) {lambda (x) (+ x (cadr pm))} 0))
(for ([pm (in-list B)])
(hash-update! h (car pm) {lambda (x) (- x (cadr pm))} 0))
(define hN (make-hash))
(define hB (make-hash))
(for ([(p m) (in-hash h)])
(cond
[(> m 0) (hash-update! hN p {lambda (x) (+ x m)} 0)]
[(< m 0) (hash-update! hB p {lambda (x) (- x m)} 0)]
[else #;(= m 0) (void)]))
(values (hash->factorization hN) (hash->factorization hB)))
(define-values (Nf-heavy B-heavy) (factorization-split Nf-factorization B-factorization))
(display "--------------------------")
(newline)
(displayln "N! heavy part excluding 2: ")
(display "Number of primes excluding 2: ")
(displayln (length Nf-heavy))
(display "Number of primes with multiplicity excluding 2: ")
(displayln (factorization-total Nf-heavy))
(display "Partial list excluding 2: ")
(displayln (append (take Nf-heavy 7)
(list "...")
(take-right Nf-heavy 3)))
(displayln "B heavy part: ")
(display "Number of primes: ")
(displayln (length Nf-heavy))
(display "Number of primes with multiplicity: ")
(displayln (factorization-total B-heavy))
(display "Partial list excluding 2: ")
(displayln (append (take B-heavy 7)
(list "...")
(take-right B-heavy 3)))
(newline)
(display "B-heavy has more primes than N!-heavy, with multiplicity excluding 2: ")
(displayln (- (factorization-total B-factorization)
(factorization-total Nf-factorization)))
(display "Should be 14891. Difference: ")
(displayln (- (- (factorization-total B-heavy)
(factorization-total Nf-heavy))
14891))
(newline)
La salida es
N! heavy part excluding 2:
Number of primes excluding 2: 23334
Number of primes with multiplicity excluding 2: 76915
Partial list excluding 2: ((11 97) (17 49) (19 65) (23 85) (29 12) (31 49) (37 32) ... (299977 1) (299983 1) (299993 1))
B heavy part:
Number of primes: 23334
Number of primes with multiplicity: 91856
Partial list excluding 2: ((3 5) (5 2) (7 4) (13 3) (43 9) (47 31) (61 1) ... (101977 48) (101987 48) (101999 48))
B-heavy has more primes than N!-heavy, with multiplicity excluding 2: 14941
Should be 14891. Difference: 50
Como era de esperar, obtuve 50 primos de más, contados con multiplicidad, que el número esperado.
Balanceando los primos en orden
Ahora intentamos reemplazar cada primo B-heavy por un primo N!-heavy n multiplicado por una potencia de 2, o sólo una potencia pura de 2. No pude entender bien el método explicado por Tao, así que creé uno propio. Simplemente apareando los primos (contados con multiplicidad) en orden.
Hay más primos (contados con multiplicidad) en B que en N!, así que rellenamos la lista con 1s. A pesar de que 1 no es primo, lo agregamos a la lista de factorización para horrorizar a los matemáticos y simplificar un poco el programa. Podemos poner los 1s al principio o al final, así que hay dos resultados.
El resultado de cada función es la cantidad de 2 no utilizados, así que 0 o un número positivo es bueno; podemos distribuirlos arbitrariamente o guardarlos para un problema futuro. También obtenemos la tabla de reemplazo de los primos B-heavy con los primos N!-heavy (contados con multiplicidad, puede incluir 1) y potencias de 2. Esto hace que la verificación sea más fácil.
(display "--------------------------")
(newline)
(define(multi-cons v n l)
(if (zero? n)
l
(multi-cons v (sub1 n) (cons v l))))
(define (do-balance/inorder Nh Bh m2)
(unless (= (factorization-total Nh) (factorization-total Bh)) (error))
(let loop ([Nh Nh] [Bh Bh] [m2 m2] [rtable '()])
(cond
[(null? Bh)
(values m2 (reverse rtable))]
[else
(define Np (caar Nh))
(define Nm (cadar Nh))
(define Bp (caar Bh))
(define Bm (cadar Bh))
(define q2 (max 0 (inexact->exact (ceiling (log (/ Bp Np) 2)))))
(unless (>= (* (expt 2 q2) Np) Bp) (error))
(define new-rtable (multi-cons (list Bp (* Np (expt 2 q2)))
(min Nm Bm)
rtable))
(cond
[(> Nm Bm)
(loop (cons (list Np (- Nm Bm)) (cdr Nh))
(cdr Bh)
(- m2 (* Bm q2))
new-rtable)]
[(< Nm Bm)
(loop (cdr Nh)
(cons (list Bp (- Bm Nm)) (cdr Bh))
(- m2 (* Nm q2))
new-rtable)]
[else
(loop (cdr Nh)
(cdr Bh)
(- m2 (* Bm q2))
new-rtable)])])))
(define (balance/inorder/first Nh Bh m2)
(define dif (- (factorization-total Bh) (factorization-total Nh)))
(define eNh (if (> dif 0)
(reverse (cons (list 1 dif) (reverse Nh)))
Nh))
(define eBh (if (< dif 0)
(reverse (cons (list 1 (- dif)) (reverse Bh)))
Bh))
(do-balance/inorder eNh eBh m2))
(define (balance/inorder/last Nh Bh m2)
(define dif (- (factorization-total Bh) (factorization-total Nh)))
(define eNh (if (> dif 0)
(cons (list 1 dif) Nh)
Nh))
(define eBh (if (< dif 0)
(cons (list 1 (- dif)) (reverse Bh))
Bh))
(do-balance/inorder eNh eBh m2))
(displayln "Balance primes / match first primes: ")
(define-values (leftover2/first table/first) (balance/inorder/first Nf-heavy B-heavy Nf-2multiplicity))
(display "leftover 2 multiplicity: ")
(displayln leftover2/first)
(displayln "Compare with 6853 or 6003")
(displayln "Partial conversion table: ")
(displayln (append (take table/first 7)
(list "...")
(take-right table/first 3)))
(newline)
(displayln "Balance primes / match last primes: ")
(define-values (leftover2/last table/last) (balance/inorder/last Nf-heavy B-heavy Nf-2multiplicity))
(display "leftover 2 multiplicity: ")
(displayln leftover2/last)
(displayln "Compare with 6853 or 6003")
(displayln "Partial conversion table: ")
(displayln (append (take table/last 7)
(list "...")
(take-right table/last 3)))
(newline)
La salida es
Balance primes / match first primes:
leftover 2 multiplicity: 6108
Compare with 6853 or 6003
Partial conversion table:
((3 11) (3 11) (3 11) (3 11) (3 11) (5 11) (5 11) ... (101999 131072) (101999 131072) (101999 131072))
Balance primes / match last primes:
leftover 2 multiplicity: 6497
Compare with 6853 or 6003
Partial conversion table:
((3 4) (3 4) (3 4) (3 4) (3 4) (5 8) (5 8) ... (101999 299977) (101999 299983) (101999 299993))
Obtuvimos 6108 y 6497 números 2 de sobra, así que el segundo método es mejor. No estoy seguro de por qué, pero es interesante que el primer método reemplaza 3 por 11 y el segundo por 4, que desperdicia menos.
Podemos comparar los resultados con los 6853 números 2 obtenidos por Tao con su método. Pero tenemos 50 primos más (ccm), así que podemos compararlos con 6853-50*17=6003. Son similares, así que voy a declarar unilateralmente una reproducción exitosa.
Verificación de las tablas de apareamiento
No me imagino un error en el código anterior, pero por si acaso, verifiquemos que las tablas de apareamiento tengan los mismos primos que los de B-heavy y N!-heavy y la potencia de 2 correcta.
Sería útil verificar el siguiente método, que es sencillo.
(define (conversion-table->factorizations t)
(define hN (make-hash))
(define hB (make-hash))
(for ([x (in-list t)])
; Reverse order of N! and B. Sorry.
(define pN (cadr x))
(define pB (car x))
(hash-add-factorization! hN (factorize/cache pN))
(hash-add-factorization! hB (factorize/cache pB)))
(define m2 (hash-ref hN 2 0))
(hash-remove! hN 2)
(values m2
(hash->factorization hN)
(hash->factorization hB)))
(display "--------------------------")
(newline)
(displayln "Balance primes / match first primes: ")
(define-values (m2/first/v Nf-heavy/first/v B-heavy/first/v)
(conversion-table->factorizations table/first))
(display "Check 2 multiplicity: ")
(displayln (= m2/first/v (- Nf-2multiplicity leftover2/first)))
(display "Check N!-heavy primes excluding 2: ")
(displayln (equal? Nf-heavy/first/v Nf-heavy))
(display "Check B-heavy primes: ")
(displayln (equal? B-heavy/first/v B-heavy))
(newline)
(displayln "Balance primes / match first primes: ")
(define-values (m2/last/v Nf-heavy/last/v B-heavy/last/v)
(conversion-table->factorizations table/last))
(display "Check 2 multiplicity: ")
(displayln (= m2/last/v (- Nf-2multiplicity leftover2/last)))
(display "Check N!-heavy primes excluding 2: ")
(displayln (equal? Nf-heavy/last/v Nf-heavy))
(display "Check B-heavy primes: ")
(displayln (equal? B-heavy/last/v B-heavy))
(newline)
La salida es
Balance primes / match first primes:
Check 2 multiplicity: #t
Check N!-heavy primes excluding 2: #t
Check B-heavy primes: #t
Balance primes / match first primes:
Check 2 multiplicity: #t
Check N!-heavy primes excluding 2: #t
Check B-heavy primes: #t
Todo #t, un éxito rotundo.
Un nuevo método óptimo para realizar éste balanceo
Ahora intentaremos encontrar un método óptimo y eficiente para realizar éste aparemeinto. Primero, la parte óptima.
En lugar de reemplazar cada primo B-pesado b por un primo N!-pesado n y una potencia de 2 o una potencia de 2 no negativa pura, tal que b <= n * 2^a, iremos en la dirección opuesta. Reemplazaremos cada primo N!-pesado n por un primo B-pesado b dividido por una potencia de 2 no negativa, tal que n >= b/2^a. Todo contado con multiplicidad. Aquí, nuevamente ampliamos la lista de primos N!-heavy con 1 que no es primo, pero que encaja perfectamente en el algoritmo y evita la necesidad de añadir otro caso especial. (También podemos extender la lista de primos B-heavy con 1s si es necesario).
La idea principal es bastante simple: reemplazamos el mayor n en la lista de N!-heavy por el mayor b/2^a posible, donde b es un primo o 1 en la lista de B-heavy y a es un entero no negativo. Esto es greedy, por lo que no necesitamos retroceder, y podemos resolverlo número por número, comenzando por el mayor. Lo difícil es demostrar que este orden da el resultado óptimo, en el sentido de que utiliza la menor potencia total de 2 posible. Tuve que afear el método para demostrar que esto es óptimo. :(
Demostración:
Por razones misteriosas que se aclaran más adelante, tengan paciencia, en la lista de primos B-heavy (o 1) permitiremos primos (o 1) divididos por una potencia de 2 no negativa. Por lo tanto, la lista incluye primos, primos divididos por una potencia positiva de 2, el número 1 y 1 dividido por una potencia positiva de 2, todo con multiplicidad. El objetivo es reemplazar cada primo (o 1) en esta lista de elementos N!-heavy n, por un número en la lista de elementos B-heavy b/2^c, tal que n >= (b/2^c)/2^a, donde c y a son no negativos, y b puede ser primo o 1, y todo contado con multiplicidad si es necesario.
Agregamos 1 (contado con multiplicidad) a la lista de elementos N-heavy o B-heavy, para obtener el mismo número de elementos (contado con multiplicidad). En ambos casos. En este paso inicial, todos son enteros y ninguno es divisible por una potencia de 2, pero queda incluido en el caso general donde la lista de elementos B-heavy puede incluir fracciones binarias.
Tomamos el número mayor de la n en la lista de números N-heavy. Tenemos que analizar dos casos.
Caso 1:
Si el mayor número b en la lista de elementos B-heavy es mayor que n, entonces b es mayor que todos los números en la lista de números N-heavy. Por lo tanto, en cualquier reemplazo, tendremos que dividir b entre 2 para usarlo. Por lo tanto, reemplazamos b por b/2 en la lista de elementos con alto contenido de B y lo intentamos de nuevo usando recursividad. (Deberíamos anotar en un rincón cuántas veces dividimos un número b por 2). Después de repetir este paso varias veces, es claro que todos los elementos en la lista B-heavy van a ser menores que n, por lo que no nos quedaremos trabados en este caso para siempre.
Caso 2:
Si el mayor número b en la lista de elementos B-heavy es menor o igual que n, entonces todos los números en la lista de elementos B-heavy son menores o iguales que n.
Supongamos que elegimos una asignación arbitraria de la lista actual de números N!-heavy a la lista actual de elementos B-heavy, dividida por potencias de 2 si es necesario. En particular, nos importa el mayor de cada lista, n y b, y asumimos que n no se asigna a b/2^a.
Entonces, n se asigna a b' (no se necesita una potencia de 2 adicional) y n' se asigna a b/2^a, por lo que, en particular, n >= b/2^a (donde a es 0 o positivo).
¿Podemos intercambiar la asignación y obtener n asignado a b y n' asignado a b'/2^a?
Como b <= n, entonces asignar n a b está bien (tampoco se necesita una potencia de 2 adicional).
Como b >= b', entonces n >= b/2^a >= b'/2^a, por lo que obtenemos que n >= b'/2^a y podemos mapear n a b'/2^a.
En lugar de mapear n a b' y n' a b/2^a, también podemos considerar mapear n a b y n' a b'/2^a. Además, a veces es posible mejorar éste mapeo y usar un número menor que a; por lo tanto, después de intercambiar, usamos el mismo valor o una potencia total menor de 2.
Para encontrar la aplicación óptima, basta con considerar sólo un apareamiento que mapea el número n más grande de los N!-heavy al número b más grande de los B-heavy. Una vez que el apareamiento del número más grande se fija (con multiplicidad), podemos intentar mapear los K-1 elementos restantes en las listas con N!-heavy y B-heavy.
Después de repetir estos dos casos suficientes veces, obtenemos una lista vacía para los números N!-heavy y B-heavy, y listo. Así que con este método greedy obtenemos una de las asignaciones óptimas.
Algoritmo:
El segundo problema es implementar esto en forma eficiente. Mantenemos la lista N!-heavy ordenada de mayor a menor, por lo que es fácil encontrar el mayor n.
Ahora debemos elegir el mayor b/2^a posible con b en la lista de números B-heavy. Tenemos que partir esta lista en dos partes.
Mantenemos los números b menores que n en una lista ordenada de mayor a menor, por lo que es fácil encontrar el mayor y eliminarlo después de usarlo.
Mantenemos los números b mayores que n en un treelist ordenada por la parte no entera del logaritmo de b en base 2, es decir, key(b) = {log_2(b)} = log_2(b) - [log_2(b)]. Almacenamos el número primo original (o 1) b en el treelist (no dividido por una potencia de 2). Si hubiéramos dividido por una potencia de 2, irían en el mismo orden, guardando el número original no tenemos que dividirlos cada vez que obtenemos un número n menor.
Los treelists son árboles RRB que simulan ser una lista, o una lista implementada secretamente como un árbol RRB. Así, tenemos con orden O(log(K)) una inserción en una posición determinada, una eliminación en una posición determinada y una búsqueda en una posición determinada.
También mantenemos el trelist ordenado, de modo que podamos usar la búsqueda binaria para obtener con orden O(log^2(K)) una inserción en el orden correcto y la búsqueda del número menor más cercano en el orden. (Con un árbol simulando ser un árbol, podríamos hacer todo esto en O(log(K)), pero tenemos que usar uno de los paquetes de Racket, así que para ahorrar unos minutos de instalación, simplemente usaremos un treelist).
En el treelist, buscamos el elemento b menor que n en ese orden y que además es el más cercano en ese orden (wrapeando en 0). Usamos la parte no entera del logaritmo en base 2, así que podemos calcular justo en el momento de uso cuántas veces tendríamos que haberlo dividirlo por 2.
Así que simplemente revisamos las dos listas y comparamos lo mejor de cada una para encontrar el número b más grande, después calculamos la potencia de 2 en el momento si es necesario, en lugar de recordar o cambiar la potencia de cada número mayor cada vez.
A medida que el número mayor n en la lista N!-pesada disminuye, sacamos los números grandes no utilizados de la lista B-heavy menor y los ponemos en la lista B-heavy mayor.
Todas estas cuentas se pueden realizar en tiempo O(K log^2(K)) (o O(K log(K)) con un árbol de verdad), donde K es el número de primos en la lista N!-pesada inicial.
Código:
(require racket/treelist)
(define (treelist->factorization h)
(sort (treelist->list h) < #:key car))
(define (treelist-factorization-total N)
(define t 0)
(for([pm (in-treelist N)])
(set! t (+ t (cadr pm))))
t)
(define (key pm)
(define l (log (car pm) 2))
(- l (floor l)))
(define (treelist-find-prev/ordered l n)
(define key-n (key n))
(cond
[(treelist-empty? l)
-1]
[(< key-n (key (treelist-first l)))
-1]
[(< (key (treelist-last l)) key-n)
(sub1 (treelist-length l))]
[else
(let loop ([first 0] [last (sub1 (treelist-length l))])
(cond
[(= (- last first) 1)
first]
[else
(define mid (quotient (+ first last) 2))
(cond
[(< key-n (key (treelist-ref l mid)))
(loop first mid)]
[else
(loop mid last)])]))]))
(define (treelist-insert/ordered l n)
(treelist-insert l (add1 (treelist-find-prev/ordered l n)) n))
(define (balance/optimal Nh Bh m2)
(define dif (- (factorization-total Bh) (factorization-total Nh)))
(let loop ([rNh (reverse (cons (list 1 dif) Nh))]
[rBh-smaller (reverse Bh)]
[tBh-bigger (treelist)]
[m2 m2]
[rtable '()])
(cond
[(null? rNh)
; And we are done ...
(values m2 (reverse rtable))]
[(and (not (null? rBh-smaller)) (> (caar rBh-smaller) (caar rNh)))
; If (caar rBh-smaller) is smaller than (caar rNh) then
; move it from rBh-smaller to tBh-bigger
(loop rNh
(cdr rBh-smaller)
(treelist-insert/ordered tBh-bigger (car rBh-smaller))
m2
rtable)]
[else
(define Np (caar rNh))
(define Nm (cadar rNh))
(define-values (Bp/s Bm/s q2/s Bv/s pos/s)
(cond
[(null? rBh-smaller)
(values #f #f #f #f #f)]
[else
(define pm (car rBh-smaller))
(values (car pm) (cadr pm) 0 (car pm) #f)]))
(define-values (Bp/b Bm/b q2/b Bv/b pos/b)
(cond
[(treelist-empty? tBh-bigger)
(values #f #f #f #f #f)]
[else
(define pos (modulo (treelist-find-prev/ordered tBh-bigger (car rNh))
(treelist-length tBh-bigger)))
(define pm (treelist-ref tBh-bigger pos))
(define Bp (car pm))
(define Bm (cadr pm))
(define q2 (inexact->exact (ceiling (log (/ Bp Np) 2))))
(unless (>= q2 0) (error))
(values Bp Bm q2 (/ Bp (expt 2 q2)) pos)]))
(cond
[(and Bv/s (or (not Bv/b) (>= Bv/s Bv/b)))
; We must remove one of the smaller ones
; No changes to m2
(unless (>= Np Bv/s) (error))
(unless (>= Np Bp/s) (error))
(unless (= q2/s 0) (error))
(define new-rtable (multi-cons (list Bp/s Np)
(min Nm Bm/s)
rtable))
(cond
[(> Nm Bm/s)
(loop (cons (list Np (- Nm Bm/s)) (cdr rNh))
(cdr rBh-smaller)
tBh-bigger
m2
new-rtable)]
[(< Nm Bm/s)
(loop (cdr rNh)
(cons (list Bp/s (- Bm/s Nm)) (cdr rBh-smaller))
tBh-bigger
m2
new-rtable)]
[else
(loop (cdr rNh)
(cdr rBh-smaller)
tBh-bigger
m2
new-rtable)])]
[else
; We must use one of the bigger ones
(unless (>= Np Bv/b) (error))
(unless (<= Np Bp/b) (error))
(unless (>= (* (expt 2 q2/b) Np) Bp/b) (error))
(unless (>= q2/b 0) (error))
(define new-rtable (multi-cons (list Bp/b (* Np (expt 2 q2/b)))
(min Nm Bm/b)
rtable))
(cond
[(> Nm Bm/b)
(loop (cons (list Np (- Nm Bm/b)) (cdr rNh))
rBh-smaller
(treelist-delete tBh-bigger pos/b)
(- m2 (* Bm/b q2/b))
new-rtable)]
[(< Nm Bm/b)
(loop (cdr rNh)
rBh-smaller
(treelist-set tBh-bigger pos/b (list Bp/b (- Bm/b Nm)))
(- m2 (* Nm q2/b))
new-rtable)]
[else
(loop (cdr rNh)
rBh-smaller
(treelist-delete tBh-bigger pos/b)
(- m2 (* Bm/b q2/b))
new-rtable)])])])))
(display "--------------------------")
(newline)
(displayln "Balance primes / match optimal primes: ")
(define-values (leftover2/optimal table/optimal) (balance/optimal Nf-heavy B-heavy Nf-2multiplicity))
(display "leftover 2 multiplicity: ")
(displayln leftover2/optimal)
(displayln "Compare with 6853 or 6003")
(displayln "Partial conversion table: ")
(displayln (append (take table/optimal 7)
(list "...")
(take-right table/optimal 3)))
(displayln "Partial conversion table (sorted): ")
(define table/optimal/sorted (sort table/optimal < #:key car))
(displayln (append (take table/optimal/sorted 7)
(list "...")
(take-right table/optimal/sorted 3)))
(define-values (m2/optimal/v Nf-heavy/optimal/v B-heavy/optimal/v)
(conversion-table->factorizations table/optimal))
(display "Check 2 multiplicity: ")
(displayln (= m2/optimal/v (- Nf-2multiplicity leftover2/optimal)))
(display "Check N!-heavy primes excluding 2: ")
(displayln (equal? Nf-heavy/optimal/v Nf-heavy))
(display "Check B-heavy primes: ")
(displayln (equal? B-heavy/optimal/v B-heavy))
(newline)
Salida:
La salida es
Balance primes / match optimal primes:
leftover 2 multiplicity: 25461
Compare with 6853 or 6003
Partial conversion table:
((101999 299993) (101999 299983) (101999 299977) (101999 299969) (101999 299951) (101999 299941) (101999 299933) ... (1039 2048) (1039 2048) (1039 2048))
Partial conversion table (sorted):
((3 4) (3 4) (3 4) (3 4) (3 4) (5 8) (5 8) ... (101999 299401) (101999 299393) (101999 299389))
Check 2 multiplicity: #t
Check N!-heavy primes excluding 2: #t
Check B-heavy primes: #t
Así que ahora tenemos 25461 números 2 sobrantes, lo cual es mayor que los seis mil y pico de los otros métodos. Mostramos la tabla de conversión usando el orden en que el algoritmo los descubrió, y también la tabla de conversión en el orden habitual. También comprobamos que la tabla de conversión tenga las factorizaciones correctas.
En resumen
Escribamos todos los pasos dentro de una sola función y usémosla en algunos ejemplos.
(define (calculate-2-balance N L A)
(define Nf2-factorization (factorize-factorial N))
(define Nf-2multiplicity (if (and (not (null? Nf2-factorization))
(= (caar Nf2-factorization) 2))
(cadar Nf2-factorization)
0))
(define Nf-factorization (if (and (not (null? Nf2-factorization))
(= (caar Nf2-factorization) 2))
(cdr Nf2-factorization)
Nf2-factorization))
(define B-factorization (factorize-B-number N L (quotient N A)))
(define-values (Nf-heavy B-heavy) (factorization-split Nf-factorization B-factorization))
(define-values (leftover2/optimal table/optimal) (balance/optimal Nf-heavy B-heavy Nf-2multiplicity))
leftover2/optimal)
(display "--------------------------")
(newline)
(displayln "Check a few examples (positive or 0 is good):")
(display "Check N=300,000, L=90,000, A=50: ")
(displayln (calculate-2-balance 300000 90000 50))
(display "Check N=300,000, L=90,000, A=150: ")
(displayln (calculate-2-balance 300000 90000 150))
(display "Check N=300,000, L=98,000, A=50: ")
(displayln (calculate-2-balance 300000 98000 50))
(display "Check N=300,000, L=98,000, A=150: ")
(displayln (calculate-2-balance 300000 98000 150))
La salida es
Check a few examples (positive or 0 is good):
Check N=300,000, L=90,000, A=50: 25461
Check N=300,000, L=90,000, A=150: 38665
Check N=300,000, L=98,000, A=50: -8188
Check N=300,000, L=98,000, A=150: 4171
Obtuvimos de nuevo 25461 potencias de 2 para el caso que analizamos con 300K!; partiendo de 90K y 50 repeticiones. Obtenemos 38665 que son más potencias de 2 usando 150 repeticiones. Así que en los dos casos obtenemos buenas descomposiciones en factores grandes. No estoy seguro de la importancia de 50 para el problema general que Tao quiere resolver.
Después de algunos intentos manuales y semimanuales, si intentamos obtener una descomposición con factores mayores que L=98K con A=50 repeticiones, obtenemos -8188 potencias de 2 sobrantes, lo que significa que el algoritmo falló. Pero si usamos A=150, obtenemos 4171 potencias de 2 sobrantes. Así que demostramos que 300! se puede descomponer en 300K factores mayores que 98K.
Lleguemos a 100K
Ahora podemos usar fuerza bruta para encontrar una solución para L=100K. Parece que las A más grandes son mejores, pero la cantidad de potencias de 2 sobrantes sube y baja, y no tengo una buena intuición ni una razón para elegir una. Así que usemos la búsqueda por fuerza bruta; lleva bastante tiempo, y por eso el código de búsqueda está comentado y la salida está copiada en el código fuente. Solo ejecutamos algunos casos interesantes al final.
(display "--------------------------")
(newline)
(displayln "Brute force check for 100000 (positive or 0 is good):")
; Too slow, uncomment to run
#;(for ([i (in-range 1 500)])
(define v (calculate-2-balance 300000 100000 i))
(when (>= v 0)
(displayln (list i v (if (>= v 0) "***" "")))))
(displayln "It should show (positive or 0 is good):")
'((283 53 ***)
(284 200 ***)
(285 58 ***)
(288 129 ***)
(299 192 ***)
(300 384 ***)
(304 95 ***)
(388 62 ***)
(389 159 ***)
(390 257 ***)
(402 105 ***)
(407 173 ***)
(418 9 ***)
(419 118 ***)
(439 52 ***)
(440 170 ***)
(441 225 ***)
)
(display "Check N=300,000, L=100,000, A=50: ")
(displayln (calculate-2-balance 300000 100000 50))
(display "Check N=300,000, L=100,000, A=150: ")
(displayln (calculate-2-balance 300000 100000 150))
(display "Check N=300,000, L=100,000, A=283: ")
(displayln (calculate-2-balance 300000 100000 283))
(display "Check N=300,000, L=100,000, A=300: ")
(displayln (calculate-2-balance 300000 100000 300))
La salida es
--------------------------
Brute force check for 100000 (positive or 0 is good):
It should show (positive or 0 is good):
'((283 53 ***)
(284 200 ***)
(285 58 ***)
(288 129 ***)
(299 192 ***)
(300 384 ***)
(304 95 ***)
(388 62 ***)
(389 159 ***)
(390 257 ***)
(402 105 ***)
(407 173 ***)
(418 9 ***)
(419 118 ***)
(439 52 ***)
(440 170 ***)
(441 225 ***))
Check N=300,000, L=100,000, A=50: -16177
Check N=300,000, L=100,000, A=150: -4103
Check N=300,000, L=100,000, A=283: 53
Check N=300,000, L=100,000, A=300: 384
:
Con L = 100K, falla (negativo) con A = 50 y con A = 150, pero obtenemos una cantidad de 2 sobrantes positiva con A = 283, que es el más chico, y también con A = 300, que proporciona una cantidad de 2 sobrantes considerable para futuras aplicaciones.
Conclusión
Ahora encontramos un método óptimo para completar el algoritmo propuesto por Tao y lo usamos para encontrar una descomposición de 300K! en 300K factores mayores que 100K, por lo que resolvimos el desafío.
Código completo
Solo para hacer más fácil cortar y pegar.
#lang racket
(require math/number-theory)
(define N 300000)
(define secret-factorize-cache (make-hash))
(define (factorize/cache n)
(hash-ref! secret-factorize-cache n (lambda () (factorize n))))
;;;;;;;;;;;;;;;;;;;;;;;;
; Factorization of N!
;;;;;;;;;;;;;;;;;;;;;;;;
(define (hash-add-factorization! h f)
(for ([pair (in-list f)])
(define p (car pair))
(define m (cadr pair))
(hash-update! h p {lambda (x) (+ x m)} 0)))
(define (hash->factorization h)
(define l (for/list ([(p m) (in-hash h)])
(list p m)))
(sort l < #:key car))
(define (factorize-factorial N)
(define h (make-hash))
(for ([n (in-range N)])
(hash-add-factorization! h (factorize/cache (add1 n))))
(hash->factorization h))
(define (factorization-total N)
(define t 0)
(for([pm (in-list N)])
(set! t (+ t (cadr pm))))
t)
(display "N! factorization calculation (N=300000): ")
(define Nf2-factorization (time (factorize-factorial N)))
(display "Number of primes including 2: ")
(displayln (length Nf2-factorization))
(display "Number of primes with multiplicity including 2: ")
(displayln (factorization-total Nf2-factorization))
(display "Multiplicity of 2: ")
(define Nf-2multiplicity (if (and (not (null? Nf2-factorization))
(= (caar Nf2-factorization) 2))
(cadar Nf2-factorization)
0))
(displayln Nf-2multiplicity)
(display "Should be 299992. Difference: ")
(displayln (- 299992 Nf-2multiplicity))
(displayln "N! factorization excluding 2: ")
(define Nf-factorization (if (and (not (null? Nf2-factorization))
(= (caar Nf2-factorization) 2))
(cdr Nf2-factorization)
Nf2-factorization))
(display "Number of primes excluding 2: ")
(displayln (length Nf-factorization))
(display "Number of primes with multiplicity excluding 2: ")
(displayln (factorization-total Nf-factorization))
(display "Partial list excluding 2: ")
(displayln (append (take Nf-factorization 7)
(list "...")
(take-right Nf-factorization 3)))
(newline)
;;;;;;;;;;;;;;;;;;;;;;;;
; Factorization of B
;;;;;;;;;;;;;;;;;;;;;;;;
(define (factorize-B-number N L size)
(define h (make-hash))
(define from (if (odd? L) L (add1 L)))
(for ([n (in-range N)])
(hash-add-factorization! h (factorize/cache (+ from (* (remainder n size) 2)))))
(hash->factorization h))
(display "B factorization calculation (L=90000, A=50): ")
(define B-factorization (time (factorize-B-number N 90000 (quotient N 50))))
(display "Number of primes: ")
(displayln (length B-factorization))
(display "Number of primes with multiplicity: ")
(displayln (factorization-total B-factorization))
(display "Partial list: ")
(displayln (append (take B-factorization 7)
(list "...")
(take-right B-factorization 3)))
(newline)
(display "B has more primes than N!, with multiplicity excluding 2: ")
(displayln (- (factorization-total B-factorization)
(factorization-total Nf-factorization)))
(display "Should be 14891. Difference: ")
(displayln (- (- (factorization-total B-factorization)
(factorization-total Nf-factorization))
14891))
(newline)
;;;;;;;;;;;;;;;;;;;;;;;;
; Split Nf Heavy and B Heavy
;;;;;;;;;;;;;;;;;;;;;;;;
(define (factorization-split N B)
(define h (make-hash))
(for ([pm (in-list N)])
(hash-update! h (car pm) {lambda (x) (+ x (cadr pm))} 0))
(for ([pm (in-list B)])
(hash-update! h (car pm) {lambda (x) (- x (cadr pm))} 0))
(define hN (make-hash))
(define hB (make-hash))
(for ([(p m) (in-hash h)])
(cond
[(> m 0) (hash-update! hN p {lambda (x) (+ x m)} 0)]
[(< m 0) (hash-update! hB p {lambda (x) (- x m)} 0)]
[else #;(= m 0) (void)]))
(values (hash->factorization hN) (hash->factorization hB)))
(define-values (Nf-heavy B-heavy) (factorization-split Nf-factorization B-factorization))
(display "--------------------------")
(newline)
(displayln "N! heavy part excluding 2: ")
(display "Number of primes excluding 2: ")
(displayln (length Nf-heavy))
(display "Number of primes with multiplicity excluding 2: ")
(displayln (factorization-total Nf-heavy))
(display "Partial list excluding 2: ")
(displayln (append (take Nf-heavy 7)
(list "...")
(take-right Nf-heavy 3)))
(displayln "B heavy part: ")
(display "Number of primes: ")
(displayln (length Nf-heavy))
(display "Number of primes with multiplicity: ")
(displayln (factorization-total B-heavy))
(display "Partial list excluding 2: ")
(displayln (append (take B-heavy 7)
(list "...")
(take-right B-heavy 3)))
(newline)
(display "B-heavy has more primes than N!-heavy, with multiplicity excluding 2: ")
(displayln (- (factorization-total B-factorization)
(factorization-total Nf-factorization)))
(display "Should be 14891. Difference: ")
(displayln (- (- (factorization-total B-heavy)
(factorization-total Nf-heavy))
14891))
(newline)
;;;;;;;;;;;;;;;;;;;;;;;;
; Balance primes in order
;;;;;;;;;;;;;;;;;;;;;;;;
(display "--------------------------")
(newline)
(define(multi-cons v n l)
(if (zero? n)
l
(multi-cons v (sub1 n) (cons v l))))
(define (do-balance/inorder Nh Bh m2)
(unless (= (factorization-total Nh) (factorization-total Bh)) (error))
(let loop ([Nh Nh] [Bh Bh] [m2 m2] [rtable '()])
(cond
[(null? Bh)
(values m2 (reverse rtable))]
[else
(define Np (caar Nh))
(define Nm (cadar Nh))
(define Bp (caar Bh))
(define Bm (cadar Bh))
(define q2 (max 0 (inexact->exact (ceiling (log (/ Bp Np) 2)))))
(unless (>= (* (expt 2 q2) Np) Bp) (error))
(define new-rtable (multi-cons (list Bp (* Np (expt 2 q2)))
(min Nm Bm)
rtable))
(cond
[(> Nm Bm)
(loop (cons (list Np (- Nm Bm)) (cdr Nh))
(cdr Bh)
(- m2 (* Bm q2))
new-rtable)]
[(< Nm Bm)
(loop (cdr Nh)
(cons (list Bp (- Bm Nm)) (cdr Bh))
(- m2 (* Nm q2))
new-rtable)]
[else
(loop (cdr Nh)
(cdr Bh)
(- m2 (* Bm q2))
new-rtable)])])))
(define (balance/inorder/first Nh Bh m2)
(define dif (- (factorization-total Bh) (factorization-total Nh)))
(define eNh (if (> dif 0)
(reverse (cons (list 1 dif) (reverse Nh)))
Nh))
(define eBh (if (< dif 0)
(reverse (cons (list 1 (- dif)) (reverse Bh)))
Bh))
(do-balance/inorder eNh eBh m2))
(define (balance/inorder/last Nh Bh m2)
(define dif (- (factorization-total Bh) (factorization-total Nh)))
(define eNh (if (> dif 0)
(cons (list 1 dif) Nh)
Nh))
(define eBh (if (< dif 0)
(cons (list 1 (- dif)) (reverse Bh))
Bh))
(do-balance/inorder eNh eBh m2))
(displayln "Balance primes / match first primes: ")
(define-values (leftover2/first table/first) (balance/inorder/first Nf-heavy B-heavy Nf-2multiplicity))
(display "leftover 2 multiplicity: ")
(displayln leftover2/first)
(displayln "Compare with 6853 or 6003")
(displayln "Partial conversion table: ")
(displayln (append (take table/first 7)
(list "...")
(take-right table/first 3)))
(newline)
(displayln "Balance primes / match last primes: ")
(define-values (leftover2/last table/last) (balance/inorder/last Nf-heavy B-heavy Nf-2multiplicity))
(display "leftover 2 multiplicity: ")
(displayln leftover2/last)
(displayln "Compare with 6853 or 6003")
(displayln "Partial conversion table: ")
(displayln (append (take table/last 7)
(list "...")
(take-right table/last 3)))
(newline)
;;;;;;;;;;;;;;;;;;;;;;;;
; Verify conversion table
;;;;;;;;;;;;;;;;;;;;;;;;
(define (conversion-table->factorizations t)
(define hN (make-hash))
(define hB (make-hash))
(for ([x (in-list t)])
; Reverse order of N! and B. Sorry.
(define pN (cadr x))
(define pB (car x))
(hash-add-factorization! hN (factorize/cache pN))
(hash-add-factorization! hB (factorize/cache pB)))
(define m2 (hash-ref hN 2 0))
(hash-remove! hN 2)
(values m2
(hash->factorization hN)
(hash->factorization hB)))
(display "--------------------------")
(newline)
(displayln "Balance primes / match first primes: ")
(define-values (m2/fisrt/v Nf-heavy/first/v B-heavy/first/v)
(conversion-table->factorizations table/first))
(display "Check 2 multiplicity: ")
(displayln (= m2/fisrt/v (- Nf-2multiplicity leftover2/first)))
(display "Check N!-heavy primes excluding 2: ")
(displayln (equal? Nf-heavy/first/v Nf-heavy))
(display "Check B-heavy primes: ")
(displayln (equal? B-heavy/first/v B-heavy))
(newline)
(displayln "Balance primes / match first primes: ")
(define-values (m2/last/v Nf-heavy/last/v B-heavy/last/v)
(conversion-table->factorizations table/last))
(display "Check 2 multiplicity: ")
(displayln (= m2/last/v (- Nf-2multiplicity leftover2/last)))
(display "Check N!-heavy primes excluding 2: ")
(displayln (equal? Nf-heavy/last/v Nf-heavy))
(display "Check B-heavy primes: ")
(displayln (equal? B-heavy/last/v B-heavy))
(newline)
;;;;;;;;;;;;;;;;;;;;;;;;
; Optimal method
;;;;;;;;;;;;;;;;;;;;;;;;
(require racket/treelist)
(define (treelist->factorization h)
(sort (treelist->list h) < #:key car))
(define (treelist-factorization-total N)
(define t 0)
(for([pm (in-treelist N)])
(set! t (+ t (cadr pm))))
t)
(define (key pm)
(define l (log (car pm) 2))
(- l (floor l)))
(define (treelist-find-prev/ordered l n)
(define key-n (key n))
(cond
[(treelist-empty? l)
-1]
[(< key-n (key (treelist-first l)))
-1]
[(< (key (treelist-last l)) key-n)
(sub1 (treelist-length l))]
[else
(let loop ([first 0] [last (sub1 (treelist-length l))])
(cond
[(= (- last first) 1)
first]
[else
(define mid (quotient (+ first last) 2))
(cond
[(< key-n (key (treelist-ref l mid)))
(loop first mid)]
[else
(loop mid last)])]))]))
(define (treelist-insert/ordered l n)
(treelist-insert l (add1 (treelist-find-prev/ordered l n)) n))
(define (balance/optimal Nh Bh m2)
(define dif (- (factorization-total Bh) (factorization-total Nh)))
(let loop ([rNh (reverse (cons (list 1 dif) Nh))]
[rBh-smaller (reverse Bh)]
[tBh-bigger (treelist)]
[m2 m2]
[rtable '()])
(cond
[(null? rNh)
; And we are done ...
(values m2 (reverse rtable))]
[(and (not (null? rBh-smaller)) (> (caar rBh-smaller) (caar rNh)))
; If (caar rBh-smaller) is smaller than (caar rNh) then
; move it from rBh-smaller to tBh-bigger
(loop rNh
(cdr rBh-smaller)
(treelist-insert/ordered tBh-bigger (car rBh-smaller))
m2
rtable)]
[else
(define Np (caar rNh))
(define Nm (cadar rNh))
(define-values (Bp/s Bm/s q2/s Bv/s pos/s)
(cond
[(null? rBh-smaller)
(values #f #f #f #f #f)]
[else
(define pm (car rBh-smaller))
(values (car pm) (cadr pm) 0 (car pm) #f)]))
(define-values (Bp/b Bm/b q2/b Bv/b pos/b)
(cond
[(treelist-empty? tBh-bigger)
(values #f #f #f #f #f)]
[else
(define pos (modulo (treelist-find-prev/ordered tBh-bigger (car rNh))
(treelist-length tBh-bigger)))
(define pm (treelist-ref tBh-bigger pos))
(define Bp (car pm))
(define Bm (cadr pm))
(define q2 (inexact->exact (ceiling (log (/ Bp Np) 2))))
(unless (>= q2 0) (error))
(values Bp Bm q2 (/ Bp (expt 2 q2)) pos)]))
(cond
[(and Bv/s (or (not Bv/b) (>= Bv/s Bv/b)))
; We must remove one of the smaller ones
; No changes to m2
(unless (>= Np Bv/s) (error))
(unless (>= Np Bp/s) (error))
(unless (= q2/s 0) (error))
(define new-rtable (multi-cons (list Bp/s Np)
(min Nm Bm/s)
rtable))
(cond
[(> Nm Bm/s)
(loop (cons (list Np (- Nm Bm/s)) (cdr rNh))
(cdr rBh-smaller)
tBh-bigger
m2
new-rtable)]
[(< Nm Bm/s)
(loop (cdr rNh)
(cons (list Bp/s (- Bm/s Nm)) (cdr rBh-smaller))
tBh-bigger
m2
new-rtable)]
[else
(loop (cdr rNh)
(cdr rBh-smaller)
tBh-bigger
m2
new-rtable)])]
[else
; We must use one of the bigger ones
(unless (>= Np Bv/b) (error))
(unless (<= Np Bp/b) (error))
(unless (>= (* (expt 2 q2/b) Np) Bp/b) (error))
(unless (>= q2/b 0) (error))
(define new-rtable (multi-cons (list Bp/b (* Np (expt 2 q2/b)))
(min Nm Bm/b)
rtable))
(cond
[(> Nm Bm/b)
(loop (cons (list Np (- Nm Bm/b)) (cdr rNh))
rBh-smaller
(treelist-delete tBh-bigger pos/b)
(- m2 (* Bm/b q2/b))
new-rtable)]
[(< Nm Bm/b)
(loop (cdr rNh)
rBh-smaller
(treelist-set tBh-bigger pos/b (list Bp/b (- Bm/b Nm)))
(- m2 (* Nm q2/b))
new-rtable)]
[else
(loop (cdr rNh)
rBh-smaller
(treelist-delete tBh-bigger pos/b)
(- m2 (* Bm/b q2/b))
new-rtable)])])])))
(display "--------------------------")
(newline)
(displayln "Balance primes / match optimal primes: ")
(define-values (leftover2/optimal table/optimal) (balance/optimal Nf-heavy B-heavy Nf-2multiplicity))
(display "leftover 2 multiplicity: ")
(displayln leftover2/optimal)
(displayln "Compare with 6853 or 6003")
(displayln "Partial conversion table: ")
(displayln (append (take table/optimal 7)
(list "...")
(take-right table/optimal 3)))
(displayln "Partial conversion table (sorted): ")
(define table/optimal/sorted (sort table/optimal < #:key car))
(displayln (append (take table/optimal/sorted 7)
(list "...")
(take-right table/optimal/sorted 3)))
(define-values (m2/optimal/v Nf-heavy/optimal/v B-heavy/optimal/v)
(conversion-table->factorizations table/optimal))
(display "Check 2 multiplicity: ")
(displayln (= m2/optimal/v (- Nf-2multiplicity leftover2/optimal)))
(display "Check N!-heavy primes excluding 2: ")
(displayln (equal? Nf-heavy/optimal/v Nf-heavy))
(display "Check B-heavy primes: ")
(displayln (equal? B-heavy/optimal/v B-heavy))
(newline)
;;;;;;;;;;;;;;;;;;;;;;;;
; All together
;;;;;;;;;;;;;;;;;;;;;;;;
(define (calculate-2-balance N L A)
(define Nf2-factorization (factorize-factorial N))
(define Nf-2multiplicity (if (and (not (null? Nf2-factorization))
(= (caar Nf2-factorization) 2))
(cadar Nf2-factorization)
0))
(define Nf-factorization (if (and (not (null? Nf2-factorization))
(= (caar Nf2-factorization) 2))
(cdr Nf2-factorization)
Nf2-factorization))
(define B-factorization (factorize-B-number N L (quotient N A)))
(define-values (Nf-heavy B-heavy) (factorization-split Nf-factorization B-factorization))
(define-values (leftover2/optimal table/optimal) (balance/optimal Nf-heavy B-heavy Nf-2multiplicity))
leftover2/optimal)
(display "--------------------------")
(newline)
(displayln "Check a few examples (positive or 0 is good):")
(display "Check N=300,000, L=90,000, A=50: ")
(displayln (calculate-2-balance 300000 90000 50))
(display "Check N=300,000, L=90,000, A=150: ")
(displayln (calculate-2-balance 300000 90000 150))
(display "Check N=300,000, L=98,000, A=50: ")
(displayln (calculate-2-balance 300000 98000 50))
(display "Check N=300,000, L=98,000, A=150: ")
(displayln (calculate-2-balance 300000 98000 150))
;;;;;;;;;;;;;;;;;;;;;;;;
; Solutions for L=100000
;;;;;;;;;;;;;;;;;;;;;;;;
(display "--------------------------")
(newline)
(displayln "Brute force check for 100000 (positive or 0 is good):")
; Too slow, uncomment to run
#;(for ([i (in-range 1 500)])
(define v (calculate-2-balance 300000 100000 i))
(when (>= v 0)
(displayln (list i v (if (>= v 0) "***" "")))))
(displayln "It should show (positive or 0 is good):")
'((283 53 ***)
(284 200 ***)
(285 58 ***)
(288 129 ***)
(299 192 ***)
(300 384 ***)
(304 95 ***)
(388 62 ***)
(389 159 ***)
(390 257 ***)
(402 105 ***)
(407 173 ***)
(418 9 ***)
(419 118 ***)
(439 52 ***)
(440 170 ***)
(441 225 ***)
)
(display "Check N=300,000, L=100,000, A=50: ")
(displayln (calculate-2-balance 300000 100000 50))
(display "Check N=300,000, L=100,000, A=150: ")
(displayln (calculate-2-balance 300000 100000 150))
(display "Check N=300,000, L=100,000, A=283: ")
(displayln (calculate-2-balance 300000 100000 283))
(display "Check N=300,000, L=100,000, A=300: ")
(displayln (calculate-2-balance 300000 100000 300))
Fin
Si usaste Ctrl+Fin, te puede interesar ver la sección sobre el caso de los 100K.
