A few days ago, Terence Tao proposed a challenge to decompose the 300K! as the product of 300K factors larger than 100K. A smaller example is decomposing 10! as the product of 10 factor greater or equal than 3.
10! = 1 * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10
10! = 3 * 3 * 4 * 4 * 4 * 5 * 5 * 6 * 6 * 7
He was able to show a decomposition of 300K! as the product of 300K factors larger than only 90K, and proposed a method to try to use larger numbers. This is part of an attempt to prove a conjecture, so the complete problem is more general, but let´s try to solve this case.
[Edit: After posting this in HN, a comment by btilly explained that Andrew Sutherland solved this case a day after it was posted. Anyway, he was using a very different method.]
We will fix N=300K that is the case analyzed by Tao in the post. The idea by Tao is to start with a odd number B that is a product of 300K odds numbers larger than 100K, and attempt to fix the difference. One important point is that 300! is even, so the idea is to get try to use the 2 to help in this task.
He defines the B-heavy primes that are the primes that appear more time in B than in N!=300!, and the N!-heavy primes that are the primes that appear more time in N!=300! than in B. The idea is to replace each B-heavy prime with bigger a N!-heavy prime n multiplied by a power of 2, or by a pure power of 2. This will increase the factor of B and we will get the requested decomposition.It's important to use as few 2 as possible, so we don't run out of them. All of this is better explained in Tao's post, so it's worth reading it.
He was able to apply this method to get a decomposition with number that are larger than 90K bu the mapping and use of 2 was not optimal, so it was not good enough for 100K. We are going to find the optimal way to do this mapping and get a decomposition with number that are larger than 100K.
A preamble
We will use Racket. First we will try to reproduce the result with 90K.
We will fix N=300K and write a memoized version of factorize, because factorization is hard and slow, and we will reuse the same number a lot.
#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!=300K!
It's much better to factorize each number an collect the factorizations, than using the builtin factorial and try to factorize the result that may take a very long time (or fail, because factorizing big numbers is hard)..
(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 multiplicitry 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)
The output is
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))
It takes like 4 second to calculate the factorization, and we get the same number of 2 than Tao. (A lot of 2.)
Factorization of B (L=90K, A=50)
To reproduce the result, we will se the same L and A than Tao, so the number B is 50 copies of the odds numbers between 90K and 102K. The calculation consider the case where the number A does not divide N, and I decked to add more odd numbers starting from 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)
The output is
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
It's faster and it takes only half a second because the factorizations are memoized, I got 50 more primes counted with multiplicity. I was using a complicated method to calculate the odds number that are repeated, but I changed it to the simple and slow method, but the 50 extra primes don't disappear. Any insight or bug report is welcome.
N!-heavy and B-heavy
Now we compare the factorization of N! and B, remove the repeated primes (counted with multiplicity).
(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)
The output is
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
Unsurprisingly, I got again 50 more primes counted with multiplicity than the expected number
Balance primes in order
Now we try to replace each B-heavy prime with bigger a N!-heavy prime n multiplied by a power of 2, or a pure power of 2. I couldn't write the method explained by Tao, so I made my own. Just match the primes (counted with multiplicity) in order.
There are more primes (counted with multiplicity) in B than in N!, so we pad the list with 1s. In spite 1 is not a prime, we will add it to the factorization lis to annoy mathematicians and make the program slightly simpler..We can put the 1s at the beginning or at the end, so there are two results.
The result of each function is the amount of unused 2, so 0 or positive is good, we can distribute them arbitrarily or save them for a future problem :). Also we get the replacement table of he heavy primes in B with the heavy primes in N! (counted with multiplicity, may include 1) and powers of 2. This makes the verification easier.
(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)
The output is
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))
We got 6108 and 6497 unused 2, so the second method is better. I'm not sure why, but it's interesting that the first method replaces 3 with 11, and the second replaces with 4 that is ess wasteful.
We can compare the results with the 6853 obtained by Tao with his method. But we have 50 more primes (cwm) so we may compare them with 6853+50*17=6003. They are similar, an I will unilaterally proclaim a successful reproduction.
Verification of the conversion tables
I can´t imagine a mistake in the previous code, but just in case let's verify the conversion tables have the same primes than B-heavy and N!-heavy and the correct power of 2.
It will be useful to verify the next method that is not straightforward.
(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)
The output is
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
All #rt, a complete success.
A new method that is optimal to make this match
Now we will try to find an optimal and efficient method to do this match. First the optimal part.
Instead of replacing each B-heavy prime b with a N!-heavy prime n and a power of 2 or a pure non negative power of 2, such that b <= n * 2^a, we will go in the other direction. We will replace each N!-heavy prime n with a B-heavy prime b divided by a non negative power of 2, such that n >= b/2^a. Everything counted with multiplicity. Here we again are extending the list of N!-heavy primes with 1 that is not a prime, but fits perfectly in the algorithm and avoids the necessity of adding yet another special case. (Also we may extend the list of B-heavy primes with 1 if neccesary.)
The main idea is quite nice and simple, we replace the biggest n in the N!-heavy list with the biggest possible b/2^a, where b is a prime or 1 in the B-heavy list and a is a non negative integer. This is greedy, so we don't need backtracking, and we can solve it one number at a time, starting from the biggest one. The hard part is proving that this order gives the optimal result, in the sense that use the less total power of 2 as possible. I had to uglify the method to prove this is optimal. :(
Proof:
For obscure reason that will be clear later, just bear with me, in the list of B-heavy primes (or 1) we will allow primes (or 1) divided by a non negative power of 2. So the list includes primes, primes divided by a positive power of 2, the number 1, and 1 divided by a positive power of 2, everything with multiplicity. So the objective is to replace each prime (or 1) in this list of N!-heavy things n with a number in the list of B-heavy things b/2^c, such that n >= (b/2^c)/2^a, where c and a are non negative, and b may be a prime or 1, and everything is counted with multiplicity if necessary.
If we have add 1 (counted with multiplicity) to the list of N!-heavy things or B-heavy things,to get the same number of elements (counted with multiplicity) in both. In this initial step all of them are integer and none of them are divided by a power of 2, but it's included in the general cases where the list of B-heavy things may include binary fractions.
We take the biggest number in the n in the list of N!-heavy numbers. We have to analyze two cases.
Case 1:
If the biggest number b in the list of B-heavy things is greater than n, then b is bigger than all the numbers in the list of N!-heavy numbers. So in any replacement we will have to divide b by 2 to use it. So we replace b by b/2 in the list of B-heavy things, and try again using recursion. (We should track somewhere else how many times we divided a number b by 2.) After repeating this step a few times, it's clear that all the things in the B-heavy list will be smaller than n, so we will not get stuck in this step forever.
Case 2:
If the biggest number b in the list of B-heavy things is smaller or equal than n, then all the number in the list of B-heavy things is smaller or equal than n.
Let's suppose we choose an arbitrary mapping from the current list of N!-heavy numbers to the current list of B-heavy things divided by powers of 2 if necessary. In particular we care about the biggest one in each list n and b, and assume n in not mapped to b/2^a.
Then n is mapped to b' (no additional power of 2 needed here) and n´ is mapped to b/2^a so in particular n>=b/2^a (where a is 0 or positive).
Can we swap the mapping and get n mapper to b and n´ mapped to b'/2^a?
Since b <= n then mapping n to b is fine (no additional power of 2 needed here neither) .
Since b>=b' then n >= b/2^a >= b'/2^a, so we get than n>=b'/2^a and we can map n to b'/2^a.
So instead of mapping n to b' and n' to b/2^a, we can also consider the mapping n to b and n' to b'/2^a. Moreover, it may be possible to improve this mapping and use a smaller number than a, so after swamping we use the same of a smaller total power of 2.
So to find the optimal mapping, it's enough to consider only the mapping that map the biggest N!-heavy numbers n to the biggest B-heavy thing b. Once the biggest number mapping is fixed (with multiplicity), we can just try to map the remaining K-1 elements in the N!-heavy and B-heavy lists.
After repeating these 2 cases enough times, we get empty list for N!-heavy and B-heavy and we are done. So with this greedy method we get one of the optimal mappings.
Algorithm:
The second problem is to implement this efficiently. We keep the N!-heavy list ordered from bigger to smaller, so it's trivial to find the biggest one n.
Now we must choose the biggest possible b/2^a with b in the list of B-heavy numbers. We have to split that list.
We keep the number that are smaller than n in a list ordered from bigger to smaller, so it's easy to find the greatest one and remove it after it's used.
We keep the number that are bigger than than n in a treelist ordered by the non integer part of the logarithm of b in base 2, i.e key(b) = {log_2(b)} = log_2(b) - [log_2(b)] We store the origianl prime number (or 1) b in the treelist (not divided by a power of 2). If we would have divided by a power of 2, hey would go in the same order, but in this way we don't actually have to divide them each time we get a smaller n.
Treelists are RRB trees pretending to be a list, or a list secretly implemented as a RRB tree. So we get O(log(K)) insertion in a detemined position, removal in a detemined position and lookup in a detemined position position.
We also will keep the treelist ordered, so we can use binary search to get O(log^2(K)) insertion in the correct order and lookup of the closest smaller number in the order. (With a tree pretending to be a tree, we could do all of this in O(log(K), but we have to use one of the packages of Racket, so to save few minutes of installation, let's just use a treelist.)
In the treelist we search the item that is smaller than n in that order and is closer in that order (this wrap up around 0). We are using the non integer part of the logarithm in base 2, so if we can calculate when we use it how many times we would have to divide it by 2.
So we just look in both lists, and compare the best of each one them to find the biggest number b after adjusting on the spot the power of 2 if neccesary, instead of tracking it or changing the power that each bigger number has each time.
As the biggest number in the N!-heavy list decrease, we will remove unused big numbers from the smaller B-heavy list and insert them in the biggest B-heavy list.
So all the processes can be done in O(K log^2(K)) time (or O(K log(K)) with a real tree, where K is the number of primes in the initial N!-heavy list.
Code:
(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)
Output:
The output is
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
So we now got 25461 leftover 2, that is bigger than the six thousand and something of the other methods. We show the conversion table using in the order the algorithm to discover them, and also the conversion table in the usual order. We also check the conversion table has the correct factorizations.
All together
Let's write all the steps inside a single function, and use it in a few examples..
(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))
The output is
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
So we got again 25461 powers of 2 for the case we analyzed with 300K!; starting from 90K and 50 repetitions. We get 38665 that are more powers of 2 using 150 repetitions instead. So in both cases we get good decompositions into large factors. I'm not sure how important is 50 for the general problem that Tao wants to solve.
After some manual and semi-manual tries, if we try to get a decomposition with factors greater than L=98K with A=50 repetitions we get -8188 leftover powers of 2, that means that the algorithm failed. But if we use A=150 we get 4171 leftover powers of 2. So we proved that 300! can be decomposed in 300K factors larger than 98K.
Let's reach 100K
Now we can use brute force to find a solution for L=100K. It looks like bigger A are better, but the number of leftover power of 2 goes up and down, and I don't have a good intuition or reason to pick one. So let's use brute force search, it takes some time so the search code is commented, and the output is copied to the source code. We only run a few interesting one at the end.
(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))
The output is
--------------------------
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
:
With L=100K it fails (negative) with A=50 and A=150, but we get positive leftover 2 with A= 283 that is the smallest one and with A=300 that gives a lot of spare 2 for future applications.
Concussion
Now found an optimal method to complete the algorithm proposed by Tao and we used it to find a decomposition of 300K! in 300K factors larger than 100K, so we solved the challenge.
Complete code
Just to make cut and paste easier.
#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))
End
In case you used Ctrl+End, you may be interested to see the section about the 100K case.
