問題141 Scheme での実装例

Scheme で書いてみた。ロジックは先日書いた通りで、単にプログラムを起こしただけです。

(define (problem-141 limit)
  (let ((lim (expt limit (/ 1 4))))
    (define (iter-k q p k sum)
      (let ((nn (* (+ (* k p p p) q) k q)))
        (if (>= nn limit)
            sum
            (iter-k q p (+ k 1) (+ sum (if (integer? (sqrt nn)) nn 0))))))
    (define (iter-p q p sum)
      (if (>= (* (+ (* p p p) q) q) limit)
          sum
          (iter-p q (+ p 1) (if (= (gcd p q) 1) (+ sum (iter-k q p 1 0)) sum))))
    (define (iter-q q sum)
      (if (>= q lim)
          sum
          (iter-q (+ q 1) (+ sum (iter-p q (+ q 1) 0)))))
    (iter-q 1 0)))

実行してみます。

gosh> (time (problem-141 1e12))
;(time (problem-141 1.0e12))
; real   1.835
; user   1.812
; sys    0.015
878454337159

問題141を解く(その2)

解法のひとつです。
幾何数列に着目するってもの考えたんだけどなあ。
いまいち無理でしたｗ

幾何数列「a, ca, c^2 a」における比率を「c」とします。
gcd(p, q) == 1 の場合、c は p/q と書くことができます。
したがって、n^2 == c^3*a^2 + a がわかります。

'a' は整数 k に対して k*q^2 の形式でなければならないため、
n^2 == (k * p^3 + q) * k * q であることがわかります。
したがって、p>q および gcd(p,q) == 1 である 3 つの整数 p, q, k を探す必要があります。
制限は q < LIMIT^(1/4) であり、p は (p^ 3+q)*q < LIMIT。

まず、この最初の部分。

幾何数列「a, ca, c^2 a」における比率を「c」とします。
gcd(p, q) == 1 の場合、c は p/q と書くことができます。
したがって、n^2 == c^3*a^2 + a がわかります。

幾何数列のどれが d, q, r になるかについては、
n^2 = d * q + r なので、r は必ず r < d かつ r < q となる。
また、d, q はどちらが最大の値を取るにしても、掛け算になるのでどちらでもかまわない。
つまり、幾何数列では一番小さい a の部分が r になる。

n^2 = d * q + r

は、

n^2 = c^2a * ca + a
n^2 + c^3 * a^2 + a

と書き換えることができる。ここまでは難しくないね。
次。

'a' は整数 k に対して k*q^2 の形式でなければならないため、
n^2 == (k * p^3 + q) * k * q であることがわかります。

ここの a を整数 k に対して置き換えるってのが、なかなか思いつけないんだよね。ここが難しい。

幾何数列「a, ca, c^2 a」における比率を「c」とします。
gcd(p, q) == 1 の場合、c は p/q と書くことができます。

これから幾何数列

a, ca, c^2a

は

a, p/q * a, p^2/q^2 * a

と書ける。
各項の公比 p/q の分母部分を消したいのと、
幾何数列だからすべての項に同じものを乗算しても公比が崩れなければ問題ないので、
a = k * q^2 って事にする。
(公比が分数だと計算がすごく面倒になるし、誤差の問題がつきまとうので、分母は消してしまいたい。)

つまり幾何数列

a, p/q * a, p^2/q^2 * a

は整数 k を用いて a = k * q ^ 2 と書けるため、

k * q ^ 2, k * p * q, k * p^2

と書き替えることができる。
さらに

n^2 = c^3 * a^2 + a

は c = p/q, a = k * q^2 から

n^2 = (k * p^3 + q) * k * q

が導かれる。
ってことでいいんかな？

最後の部分

したがって、p>q および gcd(p,q) == 1 である 3 つの整数 p, q, k を探す必要があります。
制限は q < LIMIT^(1/4) であり、p は (p^ 3+q)*q < LIMIT。

1行目はいいとして、q と p をどこまで走査すればいいかの部分は、
計算の途中で q^4 が出てきます。ただし約分するので実際に q^4 は計算しないですが、
このため LIMIT^(1/4) 未満とすれば良いのかな。たぶん。
p のほうは前述の

n^2 = (k * p^3 + q) * k * q

この式から k を省くと p は (p^ 3+q)*q < LIMIT までにすれば良さそう。
なんとなくだけど。こんな理解で良いんかな？

あとはプログラムにするだけなんですけども、これは割愛します。
完成したプログラムは爆速で解答してくれるはずです。

問題141をやってみた

メインから。
時間を計るので、ちょいゴタゴタしてます。
1兆までなのでかなり時間がかかります。
いえ、力技コードなので処理時間は気にしてないですｗ

#include <iostream>
#include <chrono>

int main()
{
    long long limit = 1000000000000;

    std::chrono::system_clock::time_point start;
    std::chrono::system_clock::time_point end;

    start = std::chrono::system_clock::now();
    long long answer = solve(limit);
    end = std::chrono::system_clock::now();
    long long elasped = std::chrono::duration_cast<std::chrono::microseconds>(end - start).count();

    std::cout << "tims: " << elasped << " answer: " << answer << std::endl;
}

次は solve。

static long long solve(long long limit)
{
    long long sum = 0;

    for (long long n = 1; n * n < limit; n++)
    {
        if (isRuishinNumber(n))
        {
            std::cout << n * n << std::endl;
            sum += n * n;
        }
    }

    return sum;
}

n * n を使うのですが n も用意しておくのは sqrt() とか使って誤差が出るのを避けるためです。

最後は isRuishinNumber。

static bool isRuishinNumber(long long n)
{
    long long nn = n * n;
    for (long long d = 2; d < n; d++)
    {
        long long q = (long long)floor((double)nn / (double)d);
        long long r = nn - d * q;

        if (r == 0)
        {
            continue;
        }

        if (d * d == q * r)
        {
            return true;
        }
    }

    return false;
}

d は 2 から n (nn の平方根) 未満まで走査します。d=1 だと r が 0 に、d=n だと d と q が n になるため累進数の条件を満たしません。
さらに、d > n を計算しても d と q が入れ替わるだけなので計算は不要です。

r は「%」を使用せず、nn - d * q と四則演算のみで求めます。こちらのほうが処理が速いです。
また、r が 0 の場合は、d, q, r は幾何数列にならないので、さっさと次の処理を実行するようにします。

d, q, r の関係は、d は n (nn の平方根)より小さいため、かならず q > d となります。
また、nn を d で割った余りが r なので、かならず d > r となります。
つまり

q > d > r > 0

という関係になります。
また、q, d, r が幾何数列にならなければいけないので、

q : d = d : r

と書くことが出来ます。公比は特に求める必要はありません。
比例式は内項の積と外項の積が等しくなります。つまり、

d * d = q * r

これが q, d, r が幾何数列になる条件となります。
ここも除算を使用すると誤差が出てしまい、正しく判定できない可能性がありますから、
整数の乗算で判定するようにしています。これならば誤差の心配はいりません。

実行結果

9
10404
16900
97344
576081
6230016
7322436
12006225
36869184
37344321
70963776
196112016
256160025
1361388609
1380568336
8534988225
9729849600
12551169024
13855173264
16394241600
123383587600
142965659664
547674002500
tims: 3681753496 answer: 878454337159

答えは 878454337159 となります。
処理時間は 1 時間ぐらいかかっていますｗ

問題に正解すると Project Euler の公式ページのフォーラムが閲覧できるようになり、みんなが解いたプログラムを拝見することができ、プログラムのロジックの解説が書いてあったり、質問したり議論もできるようになっています。計算で瞬時（数ミリ秒）で解いてしまうすごいプログラムや、独特な面白いプログラムが掲載されていたりするので、すごく勉強になります。

整数の三乗根を誤差なしで求める

三乗根は Math 標準ライブラリの Pow を使って

double val = Math.Pow((double)number, 1.0 / 3.0)

こんな風に求められるけれども、number が整数のときは
これでは誤差が出るので使い物にならない。

なのでプログラムを自分で書く！

呼び出し側。最近の .NET では main なんて書かない。

using MathLibrary;

Console.WriteLine(MathLib.CubeRoot(679L * 679L * 679L));

ライブラリなんで static ライブラリでいいじゃろ。たぶん。
本丸の処理はこんなんだ。

namespace MathLibrary
{
    public static class MathLib
    {
        private static readonly int[][] cubeRootTable =
        [
            [0, 0],
            [1, 1],
            [8, 8],
            [27, 7],
            [64, 4],
            [125, 5],
            [216, 6],
            [343, 3],
            [512, 2],
            [729, 9],
        ];

        public static long CubeRoot(long num)
        {
            List<int> result = [cubeRootTable[(int)(num % 10L)][1]];

            for (long i = num / 1000L; i > 0L; i /= 1000L)
            {
                int j;
                for (j = 1; j < 10; j++)
                {
                    if ((int)(i % 1000L) < cubeRootTable[j][0])
                    {
                        break;
                    }
                }

                result.Add(cubeRootTable[j - 1][1]);
            }

            result.Reverse();

            long val = 0L;
            foreach (int v in result)
            {
                val = 10L * val + (long)v;
            }

            // 検算
            if (val * val * val == num)
            {
                return val;
            }
            else
            {
                return 0L;
            }
        }
    }
}

3乗根が整数にならない場合は 0 を返す！（どんだけ手抜きやねんｗ)

実行してみる。

679

D:\srcs\c_sharp\MathLibrary\ConsoleApp\bin\Debug\net8.0\ConsoleApp.exe (プロセス 16328) は、コード 0 で終了しました。
デバッグが停止したときに自動的にコンソールを閉じるには、[ツール] -> [オプション] -> [デバッグ] -> [デバッグの停止時に自 動的にコンソールを閉じる] を有効にします。
このウィンドウを閉じるには、任意のキーを押してください...

誤差なんか出ない！やったね！！

約数を求める(完全版)

(define nil '())

(define (gen-primes limit)
  (let ((v (make-vector (+ limit 1) 1)))
    
    (define (set-not-prime! ini-idx)
      (define (iter i)
        (if (> i limit)
            'done
            (begin
              (vector-set! v i 0)
              (iter (+ i ini-idx)))))
      (iter (* ini-idx 2)))
    
    (define (init)
      (vector-set! v 0 0)
      (vector-set! v 1 0)
      (set-not-prime! 2))

    (define (set-prime! i)
      (if (> i (floor (sqrt limit)))
          'done
          (begin
            (set-not-prime! i)
            (set-prime! (+ i 2)))))

    (define (print-prime l h)
      (define (iter i)
        (if (> i h)
            (newline)
            (begin
              (display (vector-ref v i))
              (display " ")
              (iter (+ i 1)))))
      (iter l))

    (define (count-prime)
      (define (iter i c)
        (if (> i limit)
            c
            (iter (+ i 1) (+ c (vector-ref v i)))))
      (iter 1 0))

    (define (make-primes n)
      (let ((primes (make-vector n 0)))
        (define (iter i j)
          (cond ((>= j n) primes)
                (else
                 (let ((p (vector-ref v i)))
                   (if (= p 1)
                       (begin
                         (vector-set! primes j i)
                         (iter (+ i 1) (+ j 1)))
                       (iter (+ i 1) j ))))))
        (iter 2 0)))
    
    (begin
      (init)
      (set-prime! 3)
      (make-primes (count-prime)))))

(define (prime? n)
  (let ((primes (gen-primes 100)))
    (define (search l h)
      (let* ((m (+ l (floor (/ (- h l) 2))))
             (p (vector-ref primes m)))
        (cond ((> l h) #false)
              ((= n p) #true)
              ((< n p) (search l (- m 1)))
              (else (search (+ m 1) h)))))
    (if (< n 2)
        #false
        (search 0 (- (vector-length primes) 1)))))

(define (factorize num)
  (let ((primes (gen-primes num)))
    (define (iter n i facts)
      (if (= n 1)
          facts
          (let ((p (vector-ref primes i)))
            (if (= 0 (remainder n p))
                (iter (/ n p) i (append facts (cons p nil)))
                (iter n (+ i 1) facts)))))
    (if (prime? num)
        (list num)
        (iter num 0 nil))))

;;

(define (element-of-set? x set)
  (cond ((null? set) #f)
        ((equal? x (car set)) #t)
        (else (element-of-set? x (cdr set)))))

(define (uniq set)
  (define (iter x s)
    (if (null? x)
        s
        (iter (cdr x)
              (if (element-of-set? (car x) s)
                  s
                  (cons (car x) s)))))
  (iter set nil))

(define (subsets s)
  (if (null? s)
      (list nil)
      (let ((rest (subsets (cdr s))))
        (append rest (map (lambda (x) (cons (car s) x)) rest)))))

(define (divisor num)
  (let ((facts (factorize num)))
    (cons 1
          (sort (map (lambda (x) (apply * x)) (uniq (cdr (subsets facts))))))))

gen-prime, prime?, factorize 手続きを実装しました。
それぞれ、素数列を生成する手続き、素数判定をする手続き、素因数分解した結果をリストで返す手続きです。
素数列は大きくなるのでリストではなく、vector を使っています。また、素数の生成のアルゴリズムはエラトステネスの篩を使用しています。

divisor 手続きは、引数に数値を受け取るようにして、factorize 手続きを使うように修正します。

gosh> (divisor 12)
(1 2 3 4 6 12)
gosh> (divisor 256)
(1 2 4 8 16 32 64 128 256)
gosh> (divisor 111)
(1 3 37 111)

このように約数を求めることができます。

約数を求める

(define nil '())

(define (element-of-set? x set)
  (cond ((null? set) #f)
        ((equal? x (car set)) #t)
        (else (element-of-set? x (cdr set)))))

(define (uniq set)
  (define (iter x s)
    (if (null? x)
        s
        (iter (cdr x)
              (if (element-of-set? (car x) s)
                  s
                  (cons (car x) s)))))
  (iter set nil))

(define (subsets s)
  (if (null? s)
      (list nil)
      (let ((rest (subsets (cdr s))))
        (append rest (map (lambda (x) (cons (car s) x)) rest)))))

(define (divisor facts)
  (cons 1
        (sort (map (lambda (x) (apply * x)) (uniq (cdr (subsets facts)))))))

結論から言うと、こんなプログラムになります。
12 の約数を求めてみます。12 を素因数分解すると 2 * 2 * 3 になるので、
'(2 2 3) を divisor 手続きの引数に与えて評価します。

gosh> (divisor '(2 2 3))
(1 2 3 4 6 12)

素因数から約数を求めるロジックは、素因数リストから、すべての部分集合を求めて、
それぞれの部分集合の要素を掛け算すれば求めることができます。
各手続きは、

subsets は、集合 s の部分集合をすべて求める
uniq は集合 set から重複を取り除く
element-of-set? は集合 set に x を含むか調べる

となっています。
subsets はすべての部分集合を求めます。集合 '(a b c) の場合は、

gosh> (subsets '(a b c))
(() (c) (b) (b c) (a) (a c) (a b) (a b c))

このようになります。
並びの見た目がちょっと悪いですが、これがすべての部分集合です。
空の集合 '() と、自分自身 '(a b c) も部分集合に含むことに注意してください。

素因数分解のプログラムは別の機会に気が向いたら掲載したいと思います。

売れてないのん？

だからスクエニくんさあ・・・。

PS5 独占なんて売れなくなるからやめろって言ったやんけ。いわんこっちゃないわ。マルチプラットフォームでグローバル展開すんのがデフォやねん。はじめから PC で出せや。

MOD でティファやエアリスを・・・ｸﾞﾍﾍ・・・

・・・・・

ハッ？

PS5 独占じゃ売れないからやめれっての。今やコンシューマ機より PC のほうが圧倒的に性能高いんだし、コンシューマ機のアドバンテージなんてないんだって。はっきり言おう。PS5 じゃ AAA タイトルは無理なんよ。3D アクションゲーは性能たらんくて無理なの。4K で 60FPS でないとか見向きもされんわ。性能がミドルクラスのゲーミング PC にも届いてないわ。

さらに、PS5 と PC で同時発売とかやられるとさ、性能面で PS5 に引っ張られて、微妙やねん。コントローラの問題で UI がクソになったり。キーボードで遊ぶゲームをコントローラでやるのはいいよ。それはプレイヤーの自由だからさ。キーボード＆マウスでプレイしたい人も、コントローラーでプレイしたい人も居る。選択の自由があるのは素晴らしい事だよ。しかし PS5 コントローラ前提って話になると、その制約のせいで自由を奪われてしまうんよ。

しかも CERO とか規制厨がうるさくて、PC まで影響でるしさ。とにかく PS5 がいろいろ足を引っ張ってだな、うざいねん。PS5 大好きユーザーには悪いんだけどさ。PS5 が足枷になっとるんよ。ソニーもやる気無いし。日本市場スルーだしさ。

なので、PS5 独占で販売するのはやめなさい。わかった？スクエニくん。