extra/project-euler/508/508.factor

   1 ! Copyright (c) 2023 Cecilia Knaebchen.
   2 ! See https://factorcode.org/license.txt for BSD license.
   3 USING: kernel locals math math.functions ranges sequences
   4 project-euler.common ;
   5 IN: project-euler.508
   6
   7 ! https://projecteuler.net/problem=508
   8
   9 ! DESCRIPTION
  10 ! -----------
  11
  12 ! Consider the Gaussian integer i-1. A base i-1 representation
  13 ! of a Gaussian integer a+bi is a finite sequence of digits
  14 ! d(n-1) d(n-2) ... d(1) d(0) such that:
  15 ! a+bi = d(n-1) (i-1)^(n-1) + ... + d(1) (i-1) + d(0)
  16 ! Each d(k) is in {0,1}
  17 ! There are no leading zeros, i.e. d(n-1) != 0, unless a+bi is
  18 ! itself 0.
  19
  20 ! Here are base i-1 representations of a few Gaussian integers:
  21 ! 11+24i -> 111010110001101
  22 ! 24-11i -> 110010110011
  23 !  8+ 0i -> 111000000
  24 ! -5+ 0i -> 11001101
  25 !  0+ 0i -> 0
  26
  27 ! Remarkably, every Gaussian integer has a unique base i-1
  28 ! representation!
  29
  30 ! Define f(a+bi) as the number of 1s in the unique base i-1
  31 ! representation of a+bi. For example, f(11+24i) = 9 and
  32 ! f(24-11i) = 7.
  33
  34 ! Define B(L) as the sum of f(a+bi) for all integers a, b such
  35 ! that |a| <= L and |b| <= L. For example, B(500) = 10,795,060.
  36
  37 ! Find B(10^15) mod 1,000,000,007.
  38
  39 ! SOLUTION
  40 ! --------
  41
  42 ! f(a+bi) = sum of digits in base i-1 representation of a+bi
  43 ! Recursion for f(a+bi):
  44 !  x := (a+b) mod 2
  45 !  -> f(a+bi) = f((x-a+b)/2 + (x-a-b)/2 * i) + x
  46
  47 MEMO:: fab ( a b -- n )
  48     a b [ zero? ] both? [ 0 ] [ a b + 2 rem dup a - dup [ b + 2 / ] [ b - 2 / ] bi* fab + ] if ;
  49
  50 ! B(P,Q,R,S) := Sum(a=P..Q, b=R..S, f(a+bi))
  51 ! Recursion for B(P,Q,R,S) exists, basically four slightly
  52 ! different versions of B(-S/2,-R/2,P/2,Q/2). If summation is
  53 ! over fewer than 5000 terms, we just calculate values of f.
  54
  55 MEMO:: B ( P Q R S -- n )
  56     Q P - S R - * 5000 <
  57     [
  58         P Q [a..b] R S [a..b] [ fab ] cartesian-map [ sum ] map-sum
  59     ]
  60     [
  61         S 2 / floor neg R 2 / ceiling neg P 2 / ceiling Q 2 / floor B
  62         S 2 / floor neg R 2 / ceiling neg P 1 - 2 / ceiling Q 1 - 2 / floor 4dup [ swap - 1 + ] 2bi@ * [ B ] dip +
  63         S 1 - 2 / floor neg R 1 - 2 / ceiling neg P 2 / ceiling Q 2 / floor 4dup [ swap - 1 + ] 2bi@ * 2 * [ B ] dip +
  64         S 1 - 2 / floor neg R 1 - 2 / ceiling neg P 1 + 2 / ceiling Q 1 + 2 / floor 4dup [ swap - 1 + ] 2bi@ * [ B ] dip +
  65         + + +
  66     ] if ;
  67
  68 : euler508 ( -- answer )
  69     10 15 ^ dup [ neg ] dip 2dup B 1000000007 mod ;
  70
  71 ! [ euler508 ] time
  72 ! 19 ms
  73
  74 SOLUTION: euler508