2020年4月15日水曜日

M5 ATOM LiteとATOM Matrixを試す(開封、環境構築、表示まで)


本日到着!
技適でないバージョンが届いた方もいらっしゃるようだが、我が家はシール付きだった。
 

←例の技適シール


購入

ATOM Lite/Matrixはマイコン界隈においてESPr Developerを超える衝撃かも。
4/10にスイッチサイエンス発売ニュースを見て購入即決。
4/10の営業時間後に注文し、届いたのは本日4/14の正午くらい。
コロナ騒動で自宅待機するにはぴったりのオモチャだ。
今スイッチサイエンスの在庫を見たら品切れだったので、注文が間に合ってよかった。

今まではM5Stackは高価なので見送り、M5StickCもADCが少ないのでArduino Nano互換 + ESP8266かPi Zero Wで十分だと考えていた。
対してATOMは安価・小型かつADCも多く、上記のマイコンではカバーできない新しいWifi対応のデバイス開発の道が開けそう。

環境構築

Arduino環境で動作確認するまで多少手間取ったので、要点を以下にまとめる。
太字下線部分はハマった部分。
  • 環境
    • OS:Windows 10 64bit
    • Arduino 1.8.12 64bit (zip版)
  • USB-Serial Driver
    • CP210x_VCP_Windows.zip
    • https://m5stack.com/pages/download
      • 上記リンクの「CP2014 Driver」の横のDownload
      • スイッチサイエンスのATOMのページのクイックスタートから辿れる
      • M5公式のページではなぜかリンクされていない・・・
  • Arduino設定
    1. zipを適当な作業場所に展開
    2. Arduinoを起動し、ツールバーの「ファイル」⇒「環境設定」を開く
    3. 追加のボードマネージャのURLに以下を追加
    4. ツールバーの「スケッチ」⇒「ライブラリをインクルード」⇒「ライブラリを管理」を開く
    5. 「M5Atom」と「M5StickC」と「FastLED」を導入
      • FastLEDはMatrixでもLiteのLEDでも必須のようなので注意
    6. Serialの速度は115200bpsにする
      • 速すぎるとうまく書き込めない
    7. ボードは「ESP32 Dev Module」にする
    8. 他の設定は変更しない(周波数など)
  • Arduinoのスケッチ例
    • LEDDisplayは動かない(他の方の記事でも言及あり)
    • LEDSet/Buttom/MPU6886は動く
      • 上記はATOM Lite/Matrixの両方とも

ESP8266などでAVR以外をArduinoで開発したことがあれば、ほとんどの作業は手間取らないと思う。

以下、M5に詳しい方には常識っぽいが、ArduinoやRaspberry Piしか知らない方へのコード説明。

>#include "M5Atom.h"
 ここにM5クラスなどが入っている。

>M5.dis.drawpix(<LEDのID>, <COLOR>);
 LEDやMatrixに書き込む関数。
 LEDのIDは、Liteは0固定。Matrixは0~24の整数値で場所を指定する。
 数値と場所の対応は、USBを下側にして左上が0、右上が4、左下が20、右下が24の順。

 COLORは「GRB」それぞれが0x00~0xFFとなる。
 0x000000で全消灯、0xFFFFFFで白色の全点灯となる。
 全点灯でも目視して目が痛くはない程度。発熱も数分程度なら平気だった。
 (というより、Matrixのほうは消灯していてもUSBを繋ぐと結構温かくなる)

>M5.begin(true, false, true)
 シリアルやLCDなどの機能のON/OFFを設定する。
 M5の機能を使うためにこれを setupで呼んでおく。

>M5.update()
 ボタンなどの各種情報を更新する。

これらが使いこなせれば、あとは他のArduino基板と似たように使える。
delayのような慣れ親しんだ関数もOK。


ひとまず今日はここまで。

2020年2月27日木曜日

温室の温度をAmbientに投稿しつつ、時刻でLED照明を制御する(Ras Pi Zero⇔Arduino連動)のH/W

前回の投稿の続き。

やりたいこと
・照明の12V電源を朝7時にON、夜19時にOFFする
・温室内の4点の温度をAmbient.ioに毎分送信する

やり方
Arduino Nano
・LM61 * 4点のAnalog INを温度[℃]に変換し、1秒ごとにシリアルでPiに送る
・シリアルを監視し、「t or T」が来たらGPIOをON、「f or F」が来たらGPIOをOFFする
Raspberry Pi
・各点の温度を60秒間貯めて平均値を計算し、Ambientに送信する
・1分に1回時刻を見て、tもしくはfをシリアルで送信する

上記のような機能分担とした。
それぞれのソースコードは前回の投稿の後半を参照。近いうちに解説を書くかも。

H/Wは非常にシンプルで、次のような構成とした。
Fritzingだと綺麗に書けて気持ちよい!

青線部分は、実際にはmicroUSB-B otg + miniUSBケーブルで接続する。

実物の写真はこんな感じ。
時刻に応じて照明を消す機能はD3に割り当ててあるけど、H/Wの実装はまだ。
今はLEDでテスト中。

ブレッドボードの下側部分がLM61への配線。
D3部分は、後ほどLEDの代わりに照明への12VをON/OFFするためリレー回路を追加予定。

2020年2月25日火曜日

温室の温度をAmbientに投稿しつつ、時刻でLED照明を制御する(Ras Pi Zero⇔Arduino連動)

以前に温室の温度制御を追加したが、その温度データをリアルタイムでモニタリングしたくなった。
また、温度だけでなく日照不足も補うため、LED照明を追加した。
24時間点灯は植物の負担になるらしいので朝夕に手動でON/OFFしているが、
このON/OFFを自動化したい。(ときどき付け忘れるので…)

(↑こんな感じ 見た目はめちゃくちゃ怪しいけど、健全な草花なのでご安心を)


温度モニタリングと時間での照明ON/OFFを実現するため、
Raspberry Pi Zero WとArduino Nanoを連携させてシステムを構築した。

仕様は以下2つ。
・照明の12V電源を朝7時にON、夜19時にOFFする
・温室内の4点の温度をAmbient.ioに毎分送信する

先にAmbientの表示結果を載せる。
こんな風に4点がグラフ化できる。何か月分もさかのぼって見ることもできる優れもの。

環境検討

ESP32とかPi Zeroなら単独でもできないこともないが、
安いArduino Nano互換機を大量に持っているので2ボード構成を選択した。
わざわざ2つを繋ぐのは少しダサいけど、いくつかメリットもある。
ArduinoはAnalog INで直接読めるからLM35やLM61のような安い温度計測ICを使えるが、
Piから直接読めるセンサはお値段が高め。
ESP32はADCの誤差が大きい。また、ESP32自体がPi Zero Wより高い。
対して、Pi+Arduinoなら、Internet部分をPythonで簡単に書けるうえに安く作れる。
その代わり、PiとArduinoのSerial通信を書く必要がある。

簡単に表にまとめるとこんな感じ。
Arduino + PiPiのみESP32
主要部品Arduino Nano
Pi Zero W
LM61 *4
Pi Zero W
DS18B20 *4
ESP32
LM61 * 4
概算費用400円
1,320円
60円*4
1,960円

1,320円
320円*4
(2,540円)

2,200円
60円*4
(2,240円)
開発言語C/C++
Python
※ボード間通信有
PythonC/C++
(Nanoが互換品前提なので少しずるい気もする)

ソース

本日はソースの掲示のみ。
中身の説明とハードウェアは後日掲載します。

・Raspberry Pi
#! /usr/bin/python3
import serial as sl
import ambient
from datetime import datetime as dt

am = ambient.Ambient(<AmbientのID>, "<AmbientのKey>")

tim = dt.now()
count = 0
acc = [0.0, 0.0, 0.0, 0.0]

s = sl.Serial("/dev/ttyUSB0", 9600)
while True:
    try:
        txt = s.readline()
        ary = txt.decode().strip().split(" ")
        if len(ary) != 11:
            continue
        else:
            t = [float(x) for x in [ary[1], ary[4], ary[7], ary[10]]]
            count += 1
            for idx in range(4):
                acc[idx] += t[idx]

            if (dt.now() - tim).total_seconds() >= 60:
                t = [x / count for x in acc]
                count = 0
                acc = [0.0, 0.0, 0.0, 0.0]
                tim = dt.now()
                tmp = {'d1':t[0], 'd2':t[1], 'd3':t[2], 'd4':t[3]}
                am.send(tmp) ## Ambientへ送信

                # DC12VのON/OFF
                if 7 <= tim.hour and tim.hour <= 19:
                    s.write(b"t")
                else:
                    s.write(b"f")
    except Exception as e:
        print(e)

s.close()

・Arduino
void sendTemp();
bool serialFlag = false;
unsigned long tim = 0;
const unsigned char PWR = 3; // D3

void setup(){
  Serial.begin(9600);
  delay(20);
  pinMode(LED_BUILTIN, OUTPUT);
  pinMode(PWR, OUTPUT);
  digitalWrite(LED_BUILTIN, HIGH);
  digitalWrite(PWR, HIGH);
}

void loop(){
  if(millis() - tim > 1000){
    sendTemp();
    tim = millis();
  }
}

void serialEvent(){
  while(Serial.available()){
    char c = Serial.read();
    switch(c){
      case 'T': case 't': // AC12V ON
        digitalWrite(LED_BUILTIN, HIGH);
        digitalWrite(PWR, HIGH);
        break;
      case 'F': case 'f': // AC12V OFF
        digitalWrite(LED_BUILTIN, LOW);
        digitalWrite(PWR, LOW);
        break;
      default:
        break;
    }
  }
}

void sendTemp(){
  float t[4];
  t[0] = analogRead(A7) * 500.0/1024 - 60;
  t[1] = analogRead(A6) * 500.0/1024 - 60;
  t[2] = analogRead(A5) * 500.0/1024 - 60;
  t[3] = analogRead(A4) * 500.0/1024 - 60;
  Serial.write("S0: ");
  Serial.print(t[0]);
  Serial.write(" / S1: ");
  Serial.print(t[1]);
  Serial.write(" / S2: ");
  Serial.print(t[2]);
  Serial.write(" / S3: ");
  Serial.println(t[3]);
}

以上。続きは後編へ。

2020年2月16日日曜日

美咲フォント(misaki font)をPythonで扱う<実装>

直前の記事で美咲フォントのpng画像の扱い方を検討したが、
今度はPythonで実装を行う。
美咲フォントの詳細は直前の記事か、作成者様のサイトを参照してください。

2Byte文字1文字を入力して、美咲フォントの画像から必要な部分を切り出す。
画像を扱うのでOpenCVでやってみる。

import cv2
misaki_png_path = "[path]/misaki_mincho.png"
misaki_img = cv2.imread(misaki_png_path, 0)

def imshow(im):
    cv2.imshow("image", im)
    cv2.waitKey(0)
    cv2.destroyAllWindows()

def letter2glyph(chara):
    s = chara.encode("Shift-JIS")
    
    if s[0] < 0xA0:
        y = (s[0] - 0x81) * 2
    else:
        y = (s[0] - 0xE0 + (0x9F - 0x81 + 1)) * 2
    
    if s[1] < 0x7f:
        x = s[1] - 0x40
    elif s[1] < 0x9f:
        x = s[1] - 0x41
    else:
        x = s[1] - 0x9f
        y += 1
    
    return misaki_img[y*8:y*8+8, x*8:x*8+8]

ここで、
imshow(letter2glyph("松"))
と書くと、
「松」の画像部分を切り出して表示できる。
(8x8 pixelだからタイトルバーに比べて松の字は豆粒くらい・・・。)

簡単にコードの動作は以下の通り。
まず文字をShift-JISに変換し、s[0]=1Byte目 で s[1]=2Byte目に格納する。
そこからpng画像内の文字の位置を計算する。
(x, y)はpngファイル内の各文字の場所を指す。0始まりで、xが列、yが行。

文字コードからx,yを求める部分は、前回検討した内容が重要になる。
図の分断された4つの領域からpngの連続空間にマップするため、
x,yそれぞれに1つずつ条件分岐が必要になる。
加えて、列数を半分にする条件が1つ入る。

ここが上位1Byte:0x81~0x9F、 0xE0~0xEFの条件。
前半部分は、0始まりにするため0x81を引く。
後半部分は、0xE0を引いたうえで、前半部分の行数を足している。
最後に2倍にしているのは、列が半分になる=行が2倍になる部分を反映している。
奇数行の部分は後で。
    if s[0] < 0xA0:
        y = (s[0] - 0x81) * 2
    else:
        y = (s[0] - 0xE0 + (0x9F - 0x81 + 1)) * 2
こちらは下位1Byte:0x40~0x7E、 0x80~0xFCの条件。
0x9F未満の部分は条件そのままで、その先は次の行になる部分のため、
0x41ではなく0x9Fを引いて0列目からとし、行数も1加算して奇数行にする。
    if s[1] < 0x7F:
        x = s[1] - 0x40
    elif s[1] < 0x9F:
        x = s[1] - 0x41
    else:
        x = s[1] - 0x9F
        y += 1

文字を抜き出す部分は以上。
これで、8x8のndarrayを抜き出せる。
モノクロビットマップなので、全要素は0 or 255になっている。


おまけ

char[8]への変換と、「■」と「 」での表示。
def mat2char8(mat):
    m = (1 - mat / 255) * 2 ** np.arange(7, -1, -1)
    return np.add.reduce(m, axis=1).astype(np.uint8)
def mat2Dprint(mat):
    print("\n".join(["".join(["■" if y==0 else " " for y in x]) for x in mat]))

両方とも、実行効率改善や文字数を減らす余地がありそう。
こんな感じで遊べる。
txt = "松―ログ"
matList = [letter2glyph(c) for c in txt]
mats = np.concatenate(matList, axis=1)

mat2Dprint(mats)

charList = np.array([mat2char8(m) for m in matList])
print(charList)
表示結果↓
 ■ ■ ■                     ■ ■  
■■ ■ ■          ■■  ■■■    ■■■■ 
 ■■   ■          ■■■  ■    ■  ■ 
■■  ■   ■■■■■■■■ ■   ■   ■■   ■ 
■■  ■            ■   ■       ■  
 ■ ■ ■           ■■■■■■     ■   
 ■■■■ ■          ■        ■■    
                                
[[ 84 212  98 200 200  84 122   0]  # 松
 [  0   0   0 255   0   0   0   0]  # -
 [  0 206 114  68  68 126  64   0]  # ロ
 [ 20  30  18  98   4   8  48   0]] # グ
char[8]の結果は全角ハイフンがわかりやすい。
「グ」も、下4桁が4→8→16+32→0となっていて直感的に計算できる。

次回は、この結果をマトリクスLEDに表示して遊んでみる予定。
新幹線のインジケータを作るってわけです。

美咲フォント(misaki font)をPythonで扱う<事前検討>

(実装だけが必要の場合、こちらのページへGo!

8x8ドットの日本語ビットマップフォント、美咲フォントをPythonで取り扱ってみる。
JIS第一、第二水準をサポートで、ゴシック/明朝のTTF形式とPNG形式がある。
8x8なら、char[8]で格納でき、ダイナミック制御でセレクタを3bitにすることもでき、
データや信号線の取り回しも非常にやりやすい。
こんな素晴らしいものがフリーで公開されているとは!
公開元のウェブサイト

TTF形式も公開されているのでfonttoolsを使えば汎用的に扱えるが、
今回使うのは8x8だけなのでコンパクトなPNG形式を対象にした。  
 misaki_mincho.png ・・・ 54kB
 misaki_mincho.ttf ・・・ 1,092kB
  約20倍の差がある

このため、png画像から必要な部分を抜き出すことにする。
美咲フォントのpng画像はこんな感じの並び。
 (↑PNG版の美咲フォント明朝から一部抜粋)

この並びはShift-JISの2Byte文字部分のよう。
(参考:文字コード表 シフトJIS(Shift_JIS)より)

■Shift-JIS

Shift-JISを構成する範囲は以下の通り。
  上位1Byte:0x81~0x9F、 0xE0~0xEF
  下位1Byte:0x40~0x7E、 0x80~0xFC

字面だけではイメージしづらいので、0xFF×0xFFのドットマトリクスで示す。


何故隙間だらけかというと、ASCIIや他のコードと識別しやすくするためらしい。
特に上位1ByteがASCIIの0x20(半角スペース)や0x41(A), 0x61(a)と同じだと、
日本語なのかASCIIなのか1Byte目だけでは判別できなくなってしまうからだ。

これに対し、PNG画像の並びでいくつか注意が必要な部分がある。

➀開始点のオフセット


(※画像はhttp://charset.7jp.net/様より一部借用)

実際の文字コードは0x8140からスタートする。
それに対し、png画像は最初の8x8が対応するので、このオフセットを考慮しないといけない。

②1Byte目0xA0~0xDFと、2Byte目0x7Fの空白部分



図の通り、Shift-JISの中にはドデカい空白部分がある。
対して、美咲フォントのpng画像は空白部分を詰めている。

③png画像は行数2倍、列数1/2


2つの行頭を見比べると、pngのほうは「ぁァАА・・・亜院押」なのに対し、
Shift-JISの表は「ァА・・・院魁」となっている。
これは、pngのほうが区/点を基準とした表示なのに対し、
Shift-JISの表は16進数基準の表示になっているためである。
つまり、図で示すとおり元々の1行がpng画像では2行分になっている。


これら①②③を考慮すれば、あとは文字コード変換と簡単な画像処理でうまく扱えるはず。
続きは次回。

2020年2月15日土曜日

Project Euler(21 - 27) ※27で断念

21
In [179]: def divisors(n):
     ...:     ret = []
     ...:     for i in range(1, int(n**0.5+2), 1):
     ...:         if n%i == 0:
     ...:             if int(n/i) == i:
     ...:                 ret += [i]
     ...:             else:
     ...:                 ret += [i, int(n/i)]
     ...:     return list(set(ret))

まずは約数のユニークなリストを出す

In [180]: lst = []
     ...: for i in range(2, 10000, 1):
     ...:     ami = sum(divisors(i)) - i
     ...:     if i != ami and i == sum(divisors(ami)) - ami:
     ...:         lst += [i, ami]
     ...: int(sum(lst)/2)

Out[180]: 31626

友愛数のリストを作って合計する
重複するので2で割っておく

22
In [192]: score = 0
     ...: with open("p022_names.txt") as f:
     ...:     lst = sorted(f.readlines()[0].strip().replace('"', '').split(","))
     ...:     i = 1
     ...:     for name in lst:
     ...:         score += i * sum([ord(x)-64 for x in name])
     ...:         i+=1
     ...: score
Out[192]: 871198282

Pythonのsortedはアルファベット順に並び替えてくれるのでそのまま使える
あとは文字をforで切り出してordでASCIIコードに変えると簡潔に記述できる

23
In [197]: def isAbundant(n):
     ...:     if n < sum(divisors(n)) - n:
     ...:         return True
     ...:     else:
     ...:         return False

上記の関数でabundant number=過剰数というのを判定させる

28123 - 12以下のすべてのnに対して、過剰数リストを作る
In [236]: abundantList = []
     ...: for i in range(12, 28124, 1):
     ...:     if isAbundant(i):
     ...:         abundantList.append(i)


ついで、1から28123までで、過剰数の和で表せるものをユニークなリストにする
ここで、i = j となる12+12=24のような数も含まれるのでrangeの条件式に注意(見事にハマりました)
In [282]: lst = []
     ...: for i in range(len(abundantList)):
     ...:     for j in range(i, len(abundantList), 1):
     ...:         n = abundantList[i]+abundantList[j]
     ...:         if n < 28123:
     ...:             lst.append(n)
     ...: lst = sorted(list(set(lst)))


最後に、総数から過剰数の和で表せる数の合計を引けばOK
In [283]: sum([x for x in range(28123)]) - sum(list(set(lst)))
Out[283]: 4179871


24
先に問題を整理する
1th 0123456789 ←下1桁の入れ替えは1!=1パターン(実質入れ替えなし)
2nd 0123456798 ←下2桁の入れ替えは1th, 2ndの2!=2パターン
3rd 0123456879 
4th 0123456897
5th 0123456978
6th 0123456987 ←下3桁の入れ替えは1st-6thの3!=6パターン
つまり、n! < 1000 < (n+1)!のところを探せばかなり範囲が絞れる

このとき、n!thの数字は、下n桁が逆転した数値になる
2ndなら8,9が逆転で、3!=6thなら7,8,9が987になっている

In [291]: [factorial(x) for x in range(2,13,1)]
Out[291]: [2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600]

これより、9! - 10!の間を探せばよい
362,880thは9桁逆転なので、0987654321となる
362,881thは1023456789で、そこから8!=40,320進むとここから下8桁が逆転になる
つまり403,300thは1098765432
362,880thから9!進んだときも、同様に362,880thのうち下9桁が逆転するので、
9!*2thは1987654320となるはず

上記より、1Mthが階乗の和で表して、各桁の状態を推定する
In [313]: factorial(9)*2+factorial(8)*6+factorial(7)*6+factorial(6)*2+factorial(5)*5+factorial(4)+factorial(3)*2+factorial(2)*2
Out[313]: 1000000

つまり、9!x2+8!x6+7!x6+6!x2+5!x5+4!+3!x2+2!+1!=999,999

まず、9!x2thは1987654320, 9!x2+1th = 2013456789
ここから8!x6なので2798654310, 9!x2+8!x6+1th = 2
こうしてみると、各i!の係数aがi+1桁目の数に対応しているとわかる
n番目の数を求めるとき、n!=sigma i! x A(i)と表すと、i=10桁から順番に、i桁目には残った0-9の数字の中でA(i-1)+1番目に小さい数を当てはめていく
In [338]: n=""
     ...: nums = [0,1,2,3,4,5,6,7,8,9]
     ...: for i in [2,6,6,2,5,1,2,1,1,0]:
     ...:     n+=str(nums[i])
     ...:     nums.remove(nums[i])
     ...: int(n)
Out[338]: 2783915460

25
In [341]: a, b = 1, 1
     ...: i = 2
     ...: while len(str(b)) < 1000:
     ...:     a, b = b, a+b
     ...:     i+=1
     ...: i
Out[341]: 4782

関数を作らなくても、Whileで回すだけで十分

26
In [398]: def countDivLoop(n):
     ...:     ans = "0.0"
     ...:     lst = [10]
     ...:     divided = 100
     ...:     while True:
     ...:         a = int(divided / n)
     ...:         ans += str(a)
     ...:         divided = divided % n
     ...:         if divided == 0 or lst.count(divided) > 0:
     ...:             break
     ...:         lst.append(divided)
     ...:         divided *= 10
     ...:     return ans, lst

普通の割り算ではだめなので、割り算のループが発生するまでのリスト長さを図る関数を作成

In [403]: countDivLoop(12)
Out[403]: ('0.083', [10, 4])

In [404]: 1/12
Out[404]: 0.08333333333333333

In [405]: countDivLoop(13)
Out[405]: ('0.0769230', [10, 9, 12, 3, 4, 1])

In [406]: 1/13
Out[406]: 0.07692307692307693

In [407]: countDivLoop(14)
Out[407]: ('0.0714285', [10, 2, 6, 4, 12, 8])

In [408]: 1/14
Out[408]: 0.07142857142857142


うまく動いたので、あとは11 - 1000で最長を求める
In [412]: maxNum = 0
     ...: maxLst = []
     ...: for i in range(11, 1000, 1):
     ...:     n, lst = countDivLoop(i)
     ...:     if len(lst) > len(maxLst):
     ...:         maxNum, maxLst = i, lst
     ...: maxNum, len(maxLst)
Out[412]: (983 982)


27
In [416]: 999**2+999*999+1000
Out[416]: 1997002

この2次関数の最大値は1997002なので、それ以下の素数の列を作っておく
毎回素数判定で割り算するのはもったいないので

In [420]: def PrimeNumList(n):
     ...:     primeNumList = [2]
     ...:     i = 3
     ...:     while primeNumList[-1] < n:
     ...:         for d in primeNumList:
     ...:             if i%d == 0:
     ...:                 i+=2
     ...:                 break
     ...:             if i**0.5<d:
     ...:                 primeNumList.append(i)
     ...:                 i+=2
     ...:                 break
     ...:     return primeNumList

これに最大値を入れてリスト化しておく
これで愚直に解こうとしたけど、計算が遅すぎて断念
また練り直します

2020年1月26日日曜日

Project Euler (11 - 20) 17はスキップ

オイラーの続き


Problem 11
maxNum = 0

s="""08 02 22 97 38 15 00 40 00 75 04 05 07 78 52 12 50 77 91 08
49 49 99 40 17 81 18 57 60 87 17 40 98 43 69 48 04 56 62 00
81 49 31 73 55 79 14 29 93 71 40 67 53 88 30 03 49 13 36 65
52 70 95 23 04 60 11 42 69 24 68 56 01 32 56 71 37 02 36 91
22 31 16 71 51 67 63 89 41 92 36 54 22 40 40 28 66 33 13 80
24 47 32 60 99 03 45 02 44 75 33 53 78 36 84 20 35 17 12 50
32 98 81 28 64 23 67 10 26 38 40 67 59 54 70 66 18 38 64 70
67 26 20 68 02 62 12 20 95 63 94 39 63 08 40 91 66 49 94 21
24 55 58 05 66 73 99 26 97 17 78 78 96 83 14 88 34 89 63 72
21 36 23 09 75 00 76 44 20 45 35 14 00 61 33 97 34 31 33 95
78 17 53 28 22 75 31 67 15 94 03 80 04 62 16 14 09 53 56 92
16 39 05 42 96 35 31 47 55 58 88 24 00 17 54 24 36 29 85 57
86 56 00 48 35 71 89 07 05 44 44 37 44 60 21 58 51 54 17 58
19 80 81 68 05 94 47 69 28 73 92 13 86 52 17 77 04 89 55 40
04 52 08 83 97 35 99 16 07 97 57 32 16 26 26 79 33 27 98 66
88 36 68 87 57 62 20 72 03 46 33 67 46 55 12 32 63 93 53 69
04 42 16 73 38 25 39 11 24 94 72 18 08 46 29 32 40 62 76 36
20 69 36 41 72 30 23 88 34 62 99 69 82 67 59 85 74 04 36 16
20 73 35 29 78 31 90 01 74 31 49 71 48 86 81 16 23 57 05 54
01 70 54 71 83 51 54 69 16 92 33 48 61 43 52 01 89 19 67 48"""

mat = [[int(y) for y in x.split(" ")] for x in s.split("\n")]

def prd4(lst):
    return lst[0]*lst[1]*lst[2]*lst[3]

for r in range(20):
    for c in range(20-3):
        maxNum = max(prd4(mat[r][c:c+4]), maxNum)  # horizontal

for r in range(20-3):
    for c in range(20):
        maxNum = max(prd4([x[c] for x in mat[r:r+4]]), maxNum) # vertical

for r in range(20-3):
    for c in range(20-3):
        maxNum = max(mat[r][c]*mat[r+1][c+1]*mat[r+2][c+2]*mat[r+3][c+3], maxNum)

for r in range(20-3):
    for c in range(20-3):
        maxNum = max(mat[r+3][c]*mat[r+2][c+1]*mat[r+1][c+2]*mat[r][c+3], maxNum)

print(maxNum)

Problem 12
def divisors(n):
   ret = []
   for i in range(int(n**0.5)):
       if n%(i+1) == 0:
           ret += [i+1, int(n/(i+1))]
   return ret

i=2
triangleNum=1
while len(divisors(triangleNum)) < 500:
   triangleNum += i
   i+=1

#triangleNum = 76576500


Problem 13
和をstrにして最初の10文字をとるだけ

Problem 14
In [78]: def collatz(n):
    ...:     if n&1:
    ...:         return 3*n+1
    ...:     else:
    ...:         return int(n/2)
    ...:
    ...: def collatzList(m):
    ...:     lst = [m]
    ...:     while m != 1:
    ...:         m = collatz(m)
    ...:         lst.append(m)
    ...:     return lst
In [80]: maxLen = 0
    ...: maxNum = 0
    ...: for i in range(3, 1000000, 1):
    ...:     lst = collatzList(i)
    ...:     if maxLen < len(lst):
    ...:         maxLen = len(lst)
    ...:         maxNum = i
    ...: maxNum, maxLen
Out[80]: (837799, 525)
結構計算時間がかかるので注意

Problem 15
40P20=137846528820
20! ⇒ eval("*".join([str(x+1) for x in range(20)]))
!を定義して愚直に計算する。

Problem 16
eval("+".join(str(2**1000)))
evalを使うと記述量は削減できる。

Problem 17
英語という自然言語を扱う問題なのでパス

Problem 18
In [90]: s="""75
    ...: 95 64
    ...: 17 47 82
    ...: 18 35 87 10
    ...: 20 04 82 47 65
    ...: 19 01 23 75 03 34
    ...: 88 02 77 73 07 63 67
    ...: 99 65 04 28 06 16 70 92
    ...: 41 41 26 56 83 40 80 70 33
    ...: 41 48 72 33 47 32 37 16 94 29
    ...: 53 71 44 65 25 43 91 52 97 51 14
    ...: 70 11 33 28 77 73 17 78 39 68 17 57
    ...: 91 71 52 38 17 14 91 43 58 50 27 29 48
    ...: 63 66 04 68 89 53 67 30 73 16 69 87 40 31
    ...: 04 62 98 27 23 09 70 98 73 93 38 53 60 04 23"""
In [91]: pyramid = [[int(y) for y in x.split(" ")] for x in s.split("\n")]
In [92]: for j in range(15-1, 0, -1):
     ...:     for i in range(j):
     ...:         pyramid[j-1][i] += max(pyramid[j][i], pyramid[j][i+1])
     ...: pyramid
ピラミッドの底辺から順番に、最大値を計算していく。最小単位の3つから考えれば最大値は自明になる。

Problem 19
In [158]: n = 0
     ...: from datetime import datetime as dt
     ...: for i in range(1901, 2000+1, 1):
     ...:     for j in range(1, 12+1, 1):
     ...:         if dt(i, j, 1).weekday() == 6:
     ...:             n += 1
     ...: n
Out[158]: 171
車輪の再開発を防ぐという意味で、datetimeを使うのはギリギリ反則ではないと思う…

Problem 20
In [160]: def factorial(n):
     ...:     return eval("*".join([str(x) for x in range(n, 1, -1)]))
In [162]: eval("+".join(str(factorial(100))))
Out[162]: 648
階乗を定義しておく
とりあえず19個!

Project Euler(1 - 10)

簡単に説明すると、様々な数学的な課題にコンピュータを使ったりして解答するというパズル的なプロジェクト。
以下サイトで挑戦できる。
https://projecteuler.net/

最初の10問をPythonでやってみた。

Problem 1
In [28]: sum([x for x in range(1000) if x%5==0 or x%3==0])

Out[28]: 233168


Problem 2
In [32]: while fibList[-1] < 4000000:

    ...:     fibList.append(fibList[-1]+fibList[-2])

In [34]: sum([x for x in fibList if x&1==0])

Out[34]: 4613732


Problem 3
In [52]: n = 600851475143

    ...: lst = []

    ...: i = 3

    ...: while n >= i:

    ...:     if n%i==0:

    ...:         lst.append(i)

    ...:         n/=i

    ...:     else:

    ...:         i+=1

    ...: max(lst)

Out[52]: 6857


Problem 4
In [53]: lst = []

    ...: for i in range(100,1000,1):

    ...:     for j in range(100,1000,1):

    ...:         n = str(i*j)

    ...:         if n == "".join(list(reversed(n))):

    ...:             lst.append(int(n))

    ...: max(lst)

Out[53]: 906609


Problem 5
In [57]: 2*2*2*2*3*3*5*7*11*13*17*19

Out[57]: 232792560


Problem 6
In [59]: sum([x for x in range(101)])**2 - sum([x**2 for x in range(101)])

Out[59]: 25164150


Problem 7
In [60]: max(PrimeNumList(10001))

Out[60]: 104743

素数のリストを求めるコードは前の投稿参照。

Problem 8
In [75]: num="731671765313306249192251196744265747423553491949349698352031277450632623957831801698480186947885184385861

    ...: 56078911294949545950173795833195285320880551112540698747158523863050715693290963295227443043557668966489504452

    ...: 44523161731856403098711121722383113622298934233803081353362766142828064444866452387493035890729629049156044077

    ...: 23907138105158593079608667017242712188399879790879227492190169972088809377665727333001053367881220235421809751

    ...: 25454059475224352584907711670556013604839586446706324415722155397536978179778461740649551492908625693219784686

    ...: 22482839722413756570560574902614079729686524145351004748216637048440319989000889524345065854122758866688116427

    ...: 17147992444292823086346567481391912316282458617866458359124566529476545682848912883142607690042242190226710556

    ...: 26321111109370544217506941658960408071984038509624554443629812309878799272442849091888458015616609791913387549

    ...: 92005240636899125607176060588611646710940507754100225698315520005593572972571636269561882670428252483600823257

    ...: 530420752963450"



In [76]: maxNum = 0

    ...: for i in range(len(num) - 12):

    ...:     prd = eval("*".join(num[i:i+13]))

    ...:     maxNum = max([prd, maxNum])

    ...: maxNum

Out[76]: 23514624000


Problem 9
In [80]: for a in range(1,1000,1):

    ...:     for b in range(1,1000,1):

    ...:         c = 1000 - a - b

    ...:         if a**2 + b**2 == c**2:

    ...:             print(a, b, c)

    ...:             print(a*b*c)

    ...:             break

    ...:     else:

    ...:         continue

    ...:     break

    ...:

200 375 425

31875000


Problem 10
In [81]: s=0

    ...: for i in range(2000000):

    ...:     if isPrime(i):

    ...:         s+=i

    ...: s

Out[81]: 142913828922

2020年1月25日土曜日

Project EulerにチャレンジするためのSnippet(素数関係)

計算機で数学的な問題を解くEuler Projectという面白い学習ネタを発見。
答え合わせできるサイト: (捜索中)

まだ10個くらいしか解けていないけど、素数とか公約数公倍数とか似たようなコードを何度も使うのでSnippetを貯めていこうと思う。

手始めに、素数関係2つ。
n番目までの素数リストを返す関数と、素数判定関数。
とりあえずPythonで。

def PrimeNumList(n):
  primeNumList = [2]
  i = 3
  while len(primeNumList) < n:
    for d in primeNumList:
      if i%d == 0:
        i+=2
        break
      if i**0.5<d:
        primeNumList.append(i)
        i+=2
        break
  return primeNumList
def isPrime(n):
  if n == 2:
    return True
  elif n < 2 or n&1 == 0:
    return False
  for i in range(3, 1+int(n**0.5), 2):
    if n%i == 0:
      return False
    return True