Chronological

Posted on December 11, 2025 by rpflee

To read a WordPress blog in chronological order, simply add ?order=asc to the end of URL.

— Me@2025-12-09 01:21:29 PM

Euler problem 28.2

Posted on November 17, 2025 by rpflee

e28a :: Integer
e28a = 1 + sum [4*(n-2)^2 + 10*(n-1)
               | n <- [3, 5 .. 1001]]


e28b :: Integer
e28b = let n = 500
           sumSquares = n*(n+1)*(2*n+1) `div` 6
           sumLinear = n*(n+1) `div` 2
       in 16*sumSquares + 4*sumLinear + 4*n + 1

λ> :set +s
λ> e28a
669171001
(0.00 secs, 617,688 bytes)

λ> e28b
669171001
(0.00 secs, 82,400 bytes)
λ>

— Me@2025-11-17 12:00:22 AM

Euler problem 28.1

Posted on June 15, 2025 by rpflee

(defun e-28 ()
  (+ 1 (loop :for n :from 3 :to 1001 :by 2
             :sum (+ (* 4 (expt (- n 2) 2))
                     (* 10 (- n 1))))))

CL-USER> (e-28)
669171001

— Me@2025-06-13 11:01:45 PM

有乜嘢會令你心花怒放

Posted on May 27, 2025 by rpflee

— Me@2025-05-27 10:55:48 PM

xrandr, 2

Posted on April 20, 2025 by rpflee

xrandr -q

xrandr --output DisplayPort-3 --auto --primary --output HDMI-A-4 --auto --left-of DisplayPort-3

xrandr --output DisplayPort-3 --auto --primary --output HDMI-A-4 --off

— Me@2025-04-20 02:19:35 PM

Sapphire 16GB GDDR6

Posted on April 15, 2025 by rpflee

— Me@2025-04-15 11:52:10 AM

Euler problem 27.2.2

Posted on March 31, 2025 by rpflee

p_27 = -(2 * a - 1) * (a ^ 2 - a + 41)
  where
    n = 1000
    m = head $ filter (\x -> x ^ 2 - x + 41 > n) [1 ..]
    a = m - 1

Euler Problem 27 asks us to find a quadratic equation of the form $(n^2 + a \cdot n + b)$ that generates the most consecutive prime numbers starting from $n = 0$ , with the constraints $|a| < 1000$ and $|b| \leq 1000$ . The final answer is the product $a \cdot b$ for the winning pair. The “official” Haskell solution looks cryptic, but there’s a clever trick behind it. Let’s break it down step-by-step in a way that’s easier to follow.

A Handy Math Trick: Shifting the Quadratic

Imagine we tweak our quadratic $(n^2 + a \cdot n + b)$ by replacing $n$ with $n - x$ . When we expand it out, we get:

$(n - x)^2 + a(n - x) + b = n^2 + (a - 2x)n + (x^2 - ax + b)$

This is still a quadratic in $n$ , but now the coefficients are:

$a_{\text{new}} = a - 2x$
$b_{\text{new}} = x^2 - ax + b$

By choosing the right $x$ , we can transform any quadratic into a simpler form, like $n^2 + d$ or $n^2 + n + d$ . This is our starting point.

Which Form Works Best?

Form 1: $n^2 + d$
This one’s a dud. Plug in $n = 0, 1, 2, 3$ , and you’ll see every other number is even (e.g., $d, d+1, d+4, d+9$ ). With so many evens, it can’t produce a long run of primes—primes are mostly odd (except 2).
Form 2: $n^2 + n + d$
This is more promising. A mathematician named Rabinowitz proved that for this form to generate the longest sequence of primes, $4d - 1$ must be a special number called a Heegner number. There are only six possibilities, and after testing, $n^2 + n + 41$ comes out on top. It spits out 40 primes from $n = 0$ to $n = 39$ . The runner-up, $n^2 + n + 17$ , doesn’t do as well.

So, $n^2 + n + 41$ is our golden ticket. But there’s more to it.

The Symmetry Bonus

Here’s a cool trick with $n^2 + n + 41$ :

$(-1 - n)^2 + (-1 - n) + 41 = n^2 + n + 41$

This means if $n = 2$ gives a prime (like 47), then $n = -3$ (i.e., $-1 - 2$ ) does too (also 47). The sequence works both ways—positive and negative $n$ . Starting at $n = 0$ , we get 40 primes up to $n = 39$ , but we can also count backward with negative $n$ . Our goal is to shift this quadratic so we capture more of these primes from $n = 0$ onward, while keeping $a$ and $b$ within the problem’s limits.

Fitting It Into the Rules

Let’s start with $n^2 + n + 41$ (where $a = 1, b = 41$ ) and shift it using our trick:

New quadratic: $n^2 + (1 - 2x)n + (x^2 - x + 41)$
New $a = 1 - 2x$
New $b = x^2 - x + 41$

The problem says $|a| < 1000$ and $|b| \leq 1000$ , so:

$|1 - 2x| < 1000$
$|x^2 - x + 41| \leq 1000$

The second condition is trickier because $x^2$ grows fast. We need the biggest $x$ where $x^2 - x + 41 \leq 1000$ .

Crunching the Numbers

Let’s test some $x$ values:

$x = 31$ : $31^2 - 31 + 41 = 961 - 31 + 41 = 971$ (under 1000)
$x = 32$ : $32^2 - 32 + 41 = 1024 - 32 + 41 = 1033$ (over 1000)

So, $x = 31$ is the largest value that works. Now calculate:

$a = 1 - 2 \cdot 31 = 1 - 62 = -61$
$b = 971$ (from above)

Check the constraints:

$|a| = |-61| = 61 < 1000$ ✓
$|b| = 971 \leq 1000$ ✓

The Haskell Code Explained

That’s where this funky code comes in:

n = 1000
m = head $ filter (\x -> x ^ 2 - x + 41 > n) [1 ..]
a = m - 1

It finds $m = 32$ (first $x$ where $x^2 - x + 41 > 1000$ ). Then $x = m - 1 = 31$ .

The final line:

p_27 = -(2 * a - 1) * (a ^ 2 - a + 41)

This is a bit off. It uses $a = 31$ , but we need $a = 1 - 2x = -61$ . The correct product should be:

$a = -61$
$b = 971$
$a \cdot b = -61 \cdot 971 = -59231$

The code’s formula seems to be a typo or a leftover from a different derivation. It should just compute $(1 - 2x) \cdot (x^2 - x + 41)$ .

Why It Works

With $a = -61, b = 971$ , the quadratic $n^2 - 61n + 971$ generates 71 primes from $n = 0$ to $n = 70$ . This shift captures more of the $n^2 + n + 41$ sequence’s primes in the positive direction, thanks to the negative $a$ . The product -59231 is the answer.

Wrapping Up

This solution is sneaky—it starts with a prime-generating champ, shifts it just right, and fits the problem’s rules. It’s not the brute-force way (testing every $a$ and $b$ ), but a math shortcut that lands on the perfect answer. Pretty cool, right?

— Me@2025-03-24 12:25:20 PM

物理、工程、電腦、金融

Posted on March 8, 2025 by rpflee

— Me@2025-03-08 11:07:02 PM

Euler problem 27.1

Posted on February 23, 2025 by rpflee

(defparameter *primes-up-to-1000*
  (let ((sieve (make-array 1001
                           :element-type 'bit
                           :initial-element 0)))
    (loop for i from 2 to 1000
          when (zerop (sbit sieve i))
            do (loop for j from (* i i) to 1000 by i
                     do (setf (sbit sieve j) 1)))
    (loop for i from 2 to 1000
          when (zerop (sbit sieve i))
            collect i)))

(defun primep (n)
  (cond ((<= n 1) nil)
        ((<= n 1000) (member n *primes-up-to-1000*))
        (t (loop for p in *primes-up-to-1000*
                 while (<= (* p p) n)
                 never (zerop (mod n p))))))

(defun prime-sequence-length (a b)
  (loop for n from 0
        while (primep (+ (* n n) (* a n) b))
        count t))

(defun find-best-quadratic-coefficients
    (&optional (max-a 1000) (max-b 1000))
  (loop with best-length = 0
        with best-a = 0
        with best-b = 0
        for a from (- max-a) to max-a
        do (loop for b in *primes-up-to-1000*
                 when (> b max-b) do (loop-finish)
                   do (let ((len (prime-sequence-length
                                  a
                                  b)))
                        (when (> len best-length)
                          (setf best-length len
                                best-a a
                                best-b b))))
        finally (return (list best-length
                              best-a best-b))))

(defun euler-27 ()
  (destructuring-bind (len a b)
      (find-best-quadratic-coefficients)
    (format t
            "Sequence length: ~d, a: ~d, b: ~d~%"
            len a b)
    (* a b)))

(defun main ()
  (time (format t "Answer: ~d~%" (euler-27))))

; SLIME 2.28
CL-USER> (main)
Sequence length: 71, a: -61, b: 971
Answer: -59231
Evaluation took:
  0.103 seconds of real time
  0.103747 seconds of total run time (0.103637 user, 0.000110 system)
  100.97% CPU
  258,919,169 processor cycles
  0 bytes consed
  
NIL
CL-USER>

— Me@2025-02-23 03:48:23 PM

Euler problem 26.2.2

Posted on February 6, 2025 by rpflee

import Data.List (elemIndex, maximumBy) 
import Data.Ord (comparing)
import Data.Maybe (fromJust)

recurringCycle :: Int -> Int
recurringCycle d = remainders 1 []
  where
    remainders r rs
      | r == 0          = 0
      | s `elem` rs     = 1 + fromJust (elemIndex s rs)
      | otherwise       = remainders ((r * 10) `mod` d) (s : rs)
      where s = r `mod` d

p26 :: Int
p26 = fst $ maximumBy (comparing snd) [(n, recurringCycle n) | n <- [1,3..999], n `mod` 5 /= 0]

λ> :set +s
λ> p26
983
(0.11 secs, 29,770,480 bytes)
λ>

— Me@2025-02-06 08:29:54 AM

Euler problem 26.1

Posted on January 20, 2025 by rpflee

(defun recurring-cycle (d)
  "Calculate the length of the recurring cycle for 1/d."
  (labels ((remainders (d r rs)
             (if (eql r 0)
                 0               
                 (let ((s (mod r d)))
                   (if (member s rs)
                       (1+ (position s rs))
                       (remainders d (* 10 s) (cons s rs)))))))
    (remainders d 1 nil)))

(defun p26 ()
  "Find the number below 1000 with the longest recurring cycle in its decimal representation."
  (let* ((results (loop :for n :from 1 :to 999 
                        collect (cons n (recurring-cycle n))))
         (max-pair (reduce #'(lambda (a b)
                               (if (> (cdr a) (cdr b))
                                   a
                                   b))
                           results)))
    (car max-pair)))

; SLIME 2.28
CL-USER> (p26)
983
CL-USER>

— Me@2025-01-20 03:57:23 PM

Euler problem 25.3

Posted on January 13, 2025 by rpflee

showtime : false$
t0 : elapsed_real_time()$

t: 10^999$

j : 1$
while fib(j) < t do (
    j: j + 1
)$

j;

t1: elapsed_real_time()$

time_taken: t1-t0$

print("Time taken:", float(time_taken), "seconds");

4782

"Time taken:"0.4299999999999997"seconds"

— Me@2025-01-13 08:11:23 AM

Panjandrum

Posted on January 9, 2025 by rpflee

大架子

On the day when the Nobel Prize was announced, Louise gave me some very wise advice: “Now you have to write some unimportant papers.” I knew what she meant.

Even when scientists who receive the Nobel Prize resist the temptation to become panjandrums like Heisenberg, by issuing judgments about what others should be doing, there is a tendency for laureates to feel that one should stop doing the ordinary hard work of science, and instead go only for the next Big Thing.

I take pride in the fact that I have continued working hard on minor problems and have written a large number—literally hundreds—of unimportant papers, both before and after the Nobel Prize. Whatever other physicists may have learned from these papers, I have learned a lot. And at least I have escaped panjandrum-itis.

— Steven Weinberg: A Life in Physics

Euler problem 25.2.2

Posted on January 1, 2025 by rpflee

import System.CPUTime
import Text.Printf (printf)
import Data.List (findIndex)

time :: (Show a) => a -> IO a
time result = do
    start <- getCPUTime
    let computed = result
    end <- computed `seq` getCPUTime
    let diff = (fromIntegral (end - start)::Float)/(10^12)
    printf "Result: %s\n Time taken: %.6f seconds\n" (show computed) diff
    return computed
    
matrixMultiply :: Num a => [(a, a)] -> [(a, a)] -> [(a, a)]
matrixMultiply [(a11, a12), (a21, a22)] [(b11, b12), (b21, b22)] =
  [ (a11*b11 + a12*b21, a11*b12 + a12*b22)
  , (a21*b11 + a22*b21, a21*b12 + a22*b22) ]
 
matrixPower :: Num a => [(a, a)] -> Int -> [(a, a)]
matrixPower m 1 = m
matrixPower m n =
  let half = matrixPower m (n `div` 2)
      squared = matrixMultiply half half
  in if odd n then matrixMultiply squared m else squared
 
digitsInNumber :: (Show a, Integral a) => a -> Int
digitsInNumber = length . show
 
fibNth :: Integral a => a -> a
fibNth n 
  | n <= 2    = 1
  | otherwise = fst $ head $ matrixPower [(1, 1), (1, 0)] (fromIntegral (n - 1))
 
fibUpperBound :: Int -> Integer
fibUpperBound digitCount =
  let phi = 1.618033988749895
      logPhi = logBase 10 phi
      log5 = logBase 10 5
  in ceiling $ (fromIntegral digitCount - 1 + (log5 / 2)) / logPhi
 
binarySearchFibIndex :: Int -> Maybe Integer
binarySearchFibIndex digitCount =
  let upperBound = fibUpperBound digitCount
      binarySearch left right 
        | left > right = Nothing
        | otherwise =
            let mid = left + (right - left) `div` 2
                midFib = fibNth mid
                midDigits = digitsInNumber midFib
            in case compare midDigits digitCount of
                 EQ ->
                   let prevDigits = digitsInNumber $ fibNth (mid - 1)
                   in if prevDigits < digitCount then Just mid else binarySearch left (mid - 1)
                 LT -> binarySearch (mid + 1) right
                 GT -> binarySearch left (mid - 1)
  in binarySearch 1 upperBound

λ> time $ binarySearchFibIndex 1000
Result: Just 4782
 Time taken: 0.000852 seconds
Just 4782
λ>

λ> time $ binarySearchFibIndex 1000
Result: Just 4782
 Time taken: 0.006747 seconds
Just 4782
λ>

— Me@2024-12-25 07:00:13 AM

Euler problem 25.2.1

Posted on December 26, 2024 by rpflee

import System.CPUTime
import Text.Printf (printf)

-- | Times the execution of a pure function and prints the result.
time :: (Show a) => a -> IO a
time result = do
    start <- getCPUTime
    let computed = result
    end <- computed `seq` getCPUTime
    let diff = (fromIntegral (end - start)::Float)/(10^12)
    printf "Result: %s\n Time taken: %.6f seconds\n" (show computed) diff
    return computed

fibs :: [Integer]
fibs = 0:1:zipWith (+) fibs (tail fibs)

t :: Integer
t = 10^999
 
p25 :: Int
p25 = length w
    where
      w = takeWhile (< t) fibs

λ> time p25
Result: 4782
 Time taken: 0.000496 seconds
4782
λ>

λ> time p25
Result: 4782
 Time taken: 0.002918 seconds
4782
λ>

— Me@2024-12-25 07:00:13 AM

Euler problem 25.1.2

Posted on December 16, 2024 by rpflee

(defun matrix-multiply (a b)
  (destructuring-bind ((a11 a12) (a21 a22)) a
    (destructuring-bind ((b11 b12) (b21 b22)) b
      (list (list (+ (* a11 b11) (* a12 b21))
                  (+ (* a11 b12) (* a12 b22)))
            (list (+ (* a21 b11) (* a22 b21))
                  (+ (* a21 b12) (* a22 b22)))))))

(defun matrix-power (m n)
  (declare (type fixnum n))
  (cond ((= n 1) m)
        (t (let* ((half (matrix-power m (floor n 2)))
                  (squared (matrix-multiply half half)))
             (if (oddp n)
                 (matrix-multiply squared m)
                 squared)))))

(defun digits-in-number (n)
  "Count the number of digits in a number."
  (length (write-to-string n)))

(defun fib-nth (n)
  (declare (type fixnum n))
  (if (<= n 2)
      1
      (destructuring-bind ((result b)
                           (c d))
          (matrix-power '((1 1) (1 0)) (- n 1))
        (declare (ignore b c d))
        result)))

(defun fib-upper-bound (digit-count)
  "Calculate an upper bound for the index where Fibonacci has at least digit-count digits."
  (let* ((phi (float 1.618033988749895))
         (log-phi (log phi 10))
         (log-5 (log 5 10)))
    (ceiling (/ (+ digit-count (- 1 (/ log-5 2))) log-phi))))

(defun binary-search-fib-index (digit-count)
  "Find the index of the first Fibonacci number with the specified number of digits using binary search."
  (labels ((binary-search (left right)
             (when (<= left right)
               (let* ((mid (+ left (floor
                                    (/ (- right left) 2))))
                      (mid-fib (fib-nth mid))
                      (mid-digits (digits-in-number mid-fib)))
                 (cond ((= mid-digits digit-count) 
                        (let* ((prev-fib (fib-nth (1- mid)))
                               (prev-digits (digits-in-number
                                             prev-fib)))
                          (if (< prev-digits digit-count)
                              mid
                              (binary-search left (1- mid)))))
                       ((< mid-digits digit-count) 
                        (binary-search (1+ mid) right))
                       (t 
                        (binary-search left (1- mid))))))))
    (let* ((upper-bound (fib-upper-bound digit-count))
           (result (binary-search 1 upper-bound)))
      (or result (cons :upper-bound upper-bound)))))

; SLIME 2.28
CL-USER> (time (binary-search-fib-index 1000))
Evaluation took:
  0.001 seconds of real time
  
4782
CL-USER>

— Me@2024-12-10 07:23:47 PM

Euler problem 25.1.1

Posted on December 11, 2024 by rpflee

(defun fib (n)
  "Calculate the nth Fibonacci number using tail recursion."
  (labels ((fib-tail (a b count)
             (if (zerop count)
                 b
                 (fib-tail b (+ a b) (1- count)))))
    (if (<= n 2)
        1
        (fib-tail 1 1 (- n 2)))))

(defun matrix-multiply (a b)
  (destructuring-bind ((a11 a12) (a21 a22)) a
    (destructuring-bind ((b11 b12) (b21 b22)) b
      (list (list (+ (* a11 b11) (* a12 b21))
                  (+ (* a11 b12) (* a12 b22)))
            (list (+ (* a21 b11) (* a22 b21))
                  (+ (* a21 b12) (* a22 b22)))))))

(defun matrix-power (m n)
  (declare (type fixnum n))
  (cond ((= n 1) m)
        (t (let* ((half (matrix-power m (floor n 2)))
                  (squared (matrix-multiply half half)))
             (if (oddp n)
                 (matrix-multiply squared m)
                 squared)))))

(defun fib-nth (n)
  (declare (type fixnum n))
  (if (<= n 2)
      1
      (destructuring-bind ((result b) (c d))
          (matrix-power '((1 1) (1 0)) (- n 1))
        (declare (ignore b c d))
        result)))

(defun digits-in-number (n)
  "Count the number of digits in a number."
  (length (write-to-string n)))

(defun find-fib-index-0 (digit-count fib-f)
  "Find the index of the first Fibonacci number with the specified number of digits."
  (labels ((helper (index)
             (if (>= (digits-in-number
                      (funcall fib-f index))
                     digit-count)
                 index
                 (helper (1+ index)))))
    (helper 1)))

(defmacro find-fib-index (digit-count)
  `(find-fib-index-0 ,digit-count #'fib))
  
(defmacro find-fib-nth-index (digit-count)
  `(find-fib-index-0 ,digit-count #'fib-nth))

; SLIME 2.28
CL-USER> (time (find-fib-index 1000))
Evaluation took:
  0.456 seconds of real time
  
4782
CL-USER> (time (find-fib-nth-index 1000))
Evaluation took:
  0.080 seconds of real time
  
4782
CL-USER> (time (binary-search-fib-index 1000))
Evaluation took:
  0.001 seconds of real time
  
4782
CL-USER>

— Me@2024-12-10 07:23:47 PM

Euler problem 24.3

Posted on December 4, 2024 by rpflee

showtime : true;

digits: makelist(i,i,0,9);

p : permutations(digits)$

listify(p)[1000000];

Evaluation took 0.8800 seconds (1.0300 elapsed)
Evaluation took 0.0100 seconds (0.0100 elapsed)

[2,7,8,3,9,1,5,4,6,0]

— Me@2024-12-03 04:11:50 PM

Euler problem 24.2

Posted on November 19, 2024 by rpflee

import Data.List ( delete )

fac :: (Eq t, Num t) => t -> t
fac n = f n 1
  where
    f 0 acc = acc
    f m acc = f (m - 1) (m * acc)

perms :: Eq a => [a] -> Int -> [a] -> [a]
perms [] _ acc = reverse acc
perms xs n acc
  = perms (delete x xs) (mod n m) (x : acc)
  where 
    m = fac $ length xs - 1
    y = div n m
    x = xs !! y

p24 :: [Char]
p24 = perms "0123456789" 999999 []

λ> :set +s
λ> p24
"2783915460"
(0.02 secs, 121,968 bytes)
λ>

— Me@2024-11-18 07:22:15 PM

Euler problem 24.1.3

Posted on November 12, 2024 by rpflee

(defun fac (n)
  (labels
      ((f (m acc)
         (cond ((= m 0) acc)
               ('t (f (1- m) (* m acc))))))
    (f n 1)))

(defun string-to-list (str)
  (map 'list 'identity str))

(defun perms (ys k)
  (labels
      ((p (xs n acc)
         (if (null xs)
             (reverse acc)
             (let* ((m (fac (1- (length xs))))
                    (y (floor n m))
                    (x (nth y xs)))
               (p (delete x xs)
                  (mod n m)
                  (cons x acc))))))
    (p ys k nil)))

(concatenate 'string
             (perms (string-to-list "0123456789")
                    999999))

CL-USER> 
"2783915460"

— Me@2024-11-12 10:50:23 AM

Physics Town

continues

Category Archives: amaths-wars