カテゴリー: Math

投稿日:

pc(one liner Python Calculator)

ターミナルからPython一行コードを入力・実行

数値計算・代数計算にとってSymPyは強力な助っ人
Pythonを立ち上げることなく、ターミナルから簡単にPythonコードを一行で入力・実行できるようにするシェルスクリプト

【更新】

20230908
・インストールするファイルを1つだけにした
・モードを9つにした
・出力を7つにした
・カラー表示にした
・複数行コード入力メニュー[2c]Codeを追加

【動作確認環境】

【macOS12.6.5】
Python 3.11.4
Bash、zsh

【Ubuntu22.04.2LTS】
Python 3.11.1
Bash

【必須Pythonライブラリ】
SymPyライブラリ
mpmathライブラリ

【インストール】

HOMEにディレクトリmyscript/pcをつくる
HOME/myscript/pcに
py.sh
を配置

MacOS　.zshrcに以下を追記

source $HOME/myscript/pc/pc.sh

ubuntu　.bashrcに以下を追記

alias open=xdg-open
source $HOME/myscript/pc/pc.sh

【Menu】

[1ENTER]One-Liner [2c]Code [3f]Frac [4p]mpf [5r]Rump [6m]man [7h]history [8s]SymPyWeb [9q]Quit >
に対して1文字（ENTERなし）入力

[1][ENTER]One-Liner no change (** & ^ OK)
[2][c]Code TAB OK Enter ## to Stop inputting
[3][f]Frac Auto Change 2/3 -> Fraction(“2/3”).limit_denominator()
[4][p]mpf mp(multiple-precision) mpf Real float Auto Change 3.14 -> mpf(“3.14”)
[5][r]rump verification of Rumps example
[6][m]man This manual page
[7][h]history show code history open \$PCDIR/pclist.txt
[8][s]SymPyWeb show SymPy Web Site
[9][q]quit quit pc

【OUTPUT】

[1]pprint(eq,use_unicode=False)
[2]pprint(eq,use_unicode=True)
[3]pprint(eval(eq))
[4]print(eval(eq))
[5]N(eq, precision)
[6]latex(eval(eq))
[7]latex(N(eq,precision))

シェルスクリプト pc.sh

#!/usr/bin/env bash

# pc(one liner Python Calculator)
# Python with SymPy & Shell Scripts

version='20230908'

# 設置ディレクトリ
PCDIR=HOME/myscript/pc


# コマンド pce
# pc.sh 編集
function pce(){
openPCDIR/pc.sh
}


# コマンド pcl
# 入力履歴表示
function pcl(){
open PCDIR/pclist.txt
}

# コマンド pcd
# SymPy Doc Web表示
function pcd(){
open 'https://docs.sympy.org/latest/tutorials/intro-tutorial/features.html'
}


# コマンド help
function help(){
cat << EOF
********** General Commands Manual ***************
NAME
    pc - one liner Python Calculator

SYNTAX
    pc

VERSION
    This man page documents pc version{version}

DESCRIPTION
[Menu]
[1][ENTER]One-Liner  no change (** & ^ OK)
[2][c]Code  TAB OK Enter ## to Stop inputting
[3][f]Frac  Auto Change 2/3 -> Fraction("2/3").limit_denominator()
[4][p]mpf   mp(multiple-precision) mpf Real float  Auto Change 3.14 -> mpf("3.14")
[5][r]rump  verification of Rumps example
[6][m]man  This manual page
[7][h]history  show code history  open \PCDIR/pclist.txt
[8][s]SymPyWeb  show SymPy Web Site
[9][q]quit  quit pc

[OUTPUT] just Frac mpf
[1]pprint(eq,use_unicode=False)
[2]pprint(eq,use_unicode=True)
[3]pprint(eval(eq))
[4]print(eval(eq))
[5]N(eq, precision)
[6]latex(eval(eq))
[7]latex(N(eq,precision))

[OUTPUT] exec
[0]exec('eq')

[EXPRESSION] All SymPy Code OK
─────────────────────────┬─────────────────────────────────────────────────────
Math                     │ Python&SymPy Code
─────────────────────────┼─────────────────────────────────────────────────────
quotient                 │ //
─────────────────────────┼─────────────────────────────────────────────────────
remainder                │ %
─────────────────────────┼─────────────────────────────────────────────────────
quotient and remainder   │ divmod()
─────────────────────────┼─────────────────────────────────────────────────────
fraction                 │ 2/3 Frac(2,3) Frac("2/3")
─────────────────────────┼─────────────────────────────────────────────────────
square root              │ sqrt(2)  
─────────────────────────┼─────────────────────────────────────────────────────
Pi                       │ pi
─────────────────────────┼─────────────────────────────────────────────────────
Napier Constant          │ E
─────────────────────────┼─────────────────────────────────────────────────────
trigonometric function   │ sin(pi/3)
                         │ sin(60*pi/180)
─────────────────────────┼─────────────────────────────────────────────────────
exponential function     │ exp()
─────────────────────────┼─────────────────────────────────────────────────────
natural logarithm        │ log(E**2)
─────────────────────────┼─────────────────────────────────────────────────────
common logarithm         │ log(2,10)  log(x, base)
─────────────────────────┼─────────────────────────────────────────────────────
imaginary unit           │ (2+3j)*(5-7j)
─────────────────────────┼─────────────────────────────────────────────────────
imaginary unit(SymPy)    │ exp(cos(E**I)+sin(E*pi))
                         │ I**I
─────────────────────────┼─────────────────────────────────────────────────────
prime factorization      │ factorint(1000)
─────────────────────────┼─────────────────────────────────────────────────────
factorial                │ factorial(10)
─────────────────────────┼─────────────────────────────────────────────────────
algebra symbolic variable│ from a to z
algebra symbolic function│ only f
─────────────────────────┼─────────────────────────────────────────────────────
expand                   │ expand((x+y)**10)
─────────────────────────┼─────────────────────────────────────────────────────
factor                   │ factor(a**10-b**10)
─────────────────────────┼─────────────────────────────────────────────────────
simplification           │ simplify(1+x+x**2+x**3)
─────────────────────────┼─────────────────────────────────────────────────────
sequence                 │ sequence(k**2,(k,1,10))
                         │ sequence(k**2,(k,1,10))[9]
─────────────────────────┼─────────────────────────────────────────────────────
summation of a sequence  │ Sum(k**2,(k,1,n)).doit()
─────────────────────────┼─────────────────────────────────────────────────────
product of a sequence    │ product(k,(k,1,10))
─────────────────────────┼─────────────────────────────────────────────────────
equation                 │ solve(a*x**2+b*x+c,x)
                         │ solve(x**2-1,x)[0]
─────────────────────────┼─────────────────────────────────────────────────────
simultaneous equations   │ solve([2/3*x-y-1,3/7*x-2*y-5/9],[x,y])
─────────────────────────┼─────────────────────────────────────────────────────
simultaneous equations   │ solve([x+y-4,x-y-2],[x,y])
─────────────────────────┼─────────────────────────────────────────────────────
differential equation    │ variable function f
                         │ dsolve(Eq(f(t).diff(t, t) - f(t), exp(t)), f(t))
─────────────────────────┼─────────────────────────────────────────────────────
differential             │ diff(x**2,x)
─────────────────────────┼─────────────────────────────────────────────────────
indefinite integral      │ integrate(x**3,x)
─────────────────────────┼─────────────────────────────────────────────────────
definite integral        │ integrate(x**3,(x,0,1))
─────────────────────────┼─────────────────────────────────────────────────────
infinity                 │ oo                         │ integrate(1/(1+x**2), (x, -oo, oo))
─────────────────────────┼─────────────────────────────────────────────────────
limit                    │ limit(sin(x)/x, x, 0)
─────────────────────────┼─────────────────────────────────────────────────────
Taylor series            │ series(sin(x),x, 0, 12)
─────────────────────────┼─────────────────────────────────────────────────────
Taylor series            │ taylor(sin, 0, 5)
 coefficient list        │
─────────────────────────┼─────────────────────────────────────────────────────
substitution             │ (x - x**3/6 + x**5/120 - x**7/5040).subs(x, 1)
─────────────────────────┼─────────────────────────────────────────────────────
matrix                   │ Matrix([[1, 2], [2, 2]]).det()
─────────────────────────┼─────────────────────────────────────────────────────
matrix                   │ Matrix([[1, 2], [2, 2]]).eigenvals()
─────────────────────────┼─────────────────────────────────────────────────────
matrix                   │ Matrix([[1, 2], [2, 2]]).eigenvects()
─────────────────────────┼─────────────────────────────────────────────────────
Seki-Bernoulli number    │ bernoulli()
─────────────────────────┼─────────────────────────────────────────────────────
zeta function            │ zeta()
                         │ zetazero()
─────────────────────────┼─────────────────────────────────────────────────────
Boolean-valued check     │ 1+1 == 3                         │ expand((x+y)**2) == x**2+2*x*y+y**2
─────────────────────────┼─────────────────────────────────────────────────────
plotting graphs          │ plot(x**2, (x, -1, 2), ylabel = "y")
─────────────────────────┼─────────────────────────────────────────────────────
mpmath floating-point    │ (-2)**mpf("0.5")
─────────────────────────┴─────────────────────────────────────────────────────

********** General Commands Manual ***************

EOF
}


# コマンド pc
# 本体
function pc(){

echo -e "\033[92mpcversion\033[0m"
echo 'Copyright 2023 sakurAi Science Factory, Inc.'
echo 'This is free software with ABSOLUTELY NO WARRANTY.'


unset eq
unset N
mode=''

while :
do

while :
do
    echo -e '\033[36m[1ENTER]One-Liner [2c]Code [3f]Frac [4p]mpf [5r]Rump [6m]man [7h]history [8s]SymPyWeb [9q]Quit > \033[0m'
    stty raw -echo
    VAR=`dd bs=1 count=1 2>/dev/null`
    stty -raw echo
    case "VAR" in
    '')  mode="just" ; break ;;
    1 )  mode="just" ; break ;;
    2 )  mode="exec" ; break ;;
    c )  mode="exec" ; break ;;
    3 )  mode="frac" ; break ;;
    f )  mode="frac" ; break ;;
    4 )  mode="mpf" ; break ;;
    p )  mode="mpf" ; break ;;
    5 )  mode="rump" ; break ;;
    r )  mode="rump" ; break ;;
    6 )  help ; continue ;;
    m )  help ; continue ;;
    7 )  pcl ; continue ;;
    h )  pcl ; continue ;;
    8 )  pcd && continue ;;
    s )  pcd && continue ;;
    9 )  mode="quit" ; break ;;
    q )  mode="quit" ; break ;;
    * )  mode="just" ; break ;;
    esac
done

echo -n "">PCDIR/pceq.txt

[ "mode" = "just" ] && echo -e "\033[92m[1]eq(one liner code) >\033[0m" && read eqq && echoeqq | sed -e 's%\^%**%g' > PCDIR/pceq.txt && eq=(<PCDIR/pceq.txt) && echo -e -n '\033[92m\nprecision[default:16]>\033[0m' && read N
[ "mode" = "exec" ] && echo -e '\033[92m[2]eq(Multiple Lines + ##) >\033[0m' && eqq="" && IFS='\n' &&
while :
do
  read v
  eqq+=v'\n'
  if [ v = "##" ]; then
    break
  fi
done
echo -e "eqq" > PCDIR/pceq.txt
[ "mode" = "frac" ] && echo -e "\033[92m[3]eq(one liner code) >\033[0m" && read eqq && echo "eqq" >>PCDIR/pclist.txt && echo eqq | sed -e 's%[0-9]\+/[0-9]\+%Frac("&").limit_denominator()%g' -e 's%\^%**%g'>PCDIR/pceq.txt && eq=(<PCDIR/pceq.txt) && echo -e -n '\033[92m\nprecision[default:16] > \033[0m' && read N
[ "mode" = "mpf" ] && echo -e "\033[92m[4]eq(one liner code) >\033[0m" && read eqq && echoeqq | sed -e 's%[0-9]\+\.*[0-9]*%mpf("&")%g' -e 's%\^%**%g' > PCDIR/pceq.txt && eq=(<PCDIR/pceq.txt) && echo -e -n '\033[92m\nprecision[default:16]>\033[0m' && read N
[ "mode" = "rump" ] && echo -e '\033[92mRump'\''s example Test\033[0m' && eq="rump"
[ "mode" = "quit" ] && echo -e '\033[92mQuit pc\033[0m' && break


COMMAND=(cat << EOF
echo "eq" >>PCDIR/pclist.txt
echo -E "
from fractions import Fraction as Frac
from mpmath import *
from pprint import pprint
mp.pretty = True
from sympy import *

init_printing() #降べきの順
# init_printing(order='rev-lex') #昇べきの順
var('a:z')
f = Function('f')

GRE='\033[92m'
RED='\033[91m'
VIO='\033[95m'
END='\033[0m'

mp.dps = {N:-16}

print(VIO+'\n[1]pprint(eq,use_unicode=False)'+END)
pprint(eq,use_unicode=False)

print(VIO+'\n[2]pprint(eq,use_unicode=True)'+END)
pprint(eq,use_unicode=True)

print(VIO+'\n[3]pprint(eval(eq), order="rev-lex")'+END)
pprint(eval('eq'), order='rev-lex')

print(VIO+'\n[4]print(eval(eq))'+END)
print(eval('eq'))

print(VIO+'\n[5]N(eq,{N:-16})'+END)
try:
    N(eq,{N:- 16})
except AttributeError:
    print(RED+'No Numerical Evaluation'+END)
else:
    pprint(N(eq,{N:-16}))

print(VIO+'\n[6]latex(eval(eq))'+END)
try:
    latex(eval('eq'))
except AttributeError:
    print(RED+'No Output in LaTeX'+END)
else:
    print(latex(eval('eq')))

print(VIO+'\n[7]latex(N(eq,{N}))'+END)
try:
    latex(N(eq,{N:-16}))
except AttributeError:
    print(RED+'No Output in LaTeX'+END)
else:
    print(latex(N(eq,{N:-16})))
print('\n')
" | python

EOF

)


# モード毎
casemode in
just )
eval "{COMMAND}"
;;


exec )
echo "eqq" >> PCDIR/pclist.txt
echo -E "
from mpmath import *
from pprint import pprint
mp.pretty = True
from sympy import *
from fractions import Fraction as Frac
init_printing()
#init_printing(order='rev-lex')
var('a:z')
f = Function('f')

GRE='\033[92m'
RED='\033[91m'
VIO='\033[95m'
END='\033[0m'

mp.dps ={N:- 16}

print(VIO+'\n[0]exec(eq)'+END)
f = open('PCDIR/pceq.txt')
cmd = f.read()
exec(cmd)
print('\n')
" | python
echo -n "">PCDIR/pceq.txt
;;


frac )
echo eq && echo -n -e '\033[92mexpression change? [n(ENTER)/y]\033[0m ' ; read yn
case "yn" in [Yy])
read eq ;;
[])
;;
[n])
;;
esac
eval echo "COMMAND"
;;


mpf )
echo eq && echo -n -e '\033[92mexpression change? [n(ENTER)/y]\033[0m ' ; read yn
case "yn" in [Yy])
read eq ;;
[])
;;
[n])
;;
esac
eval "{COMMAND}"
;;


rump )
echo -e -n '\033[92mprecision[default:16]>\033[0m'
read N
echo "eq" >> PCDIR/pclist.txt
echo -E "
from mpmath import *
mp.pretty = True

a=77617
b=33096
c=333.75*b**6+a**2*(11*a**2*b**2-b**6-121*b**4-2)+5.5*b**8+a/(2*b)
print(f'Normal {c}')

def g(a, b):
    return (mpf('333.75')*b**6 + a**2*(11*a**2*b**2-b**6-121*b**4-2)+mpf('5.5')*b**8+a/(mpf('2')*b))
print('{:6}'.format('mp.dps'))
for mp.dps in range(1,N+1):
    print('{:6}'.format(mp.dps),g(mpf('77617'), mpf('33096')))
print('')
" | python

esac

done

}
投稿日:

One liner SymPy Calculator（一行数式SymPy計算機）Ver.20230508

ターミナルだけですぐに計算できる

基本は一行のPythonSymPyコードをpythonに渡すだけの仕組み
SymPyとmpmathライブラリの関数を使った一行コードが実行できる
以下仕様のようにもろもろのセッティングが不必要
いちいちpythonを起動、コーディングをしなくても計算できることを追求

【条件】

Bash、zshでの使用OK
python実行環境インストール済み
SymPyライブラリ
mpmathライブラリ

実装環境
macOS12.6.5
Python 3.9.16

Ubuntu22.04.2LTS
Python 3.8.16

【仕様】

入力数式はpython文法
シンボリック変数はaからz
シンボリック関数はf 微分方程式用
sympyとmpmathの関数はすべて使える
分数 2/3 を Frac(2,3)またはFrac(2/3)とすれば分数として計算・出力
おまけ　Rumpの例題の検証　入力にrump

【使い方】

ターミナルで

$ pycalc

で本体起動

数値計算式・代数式を入力（変数eq）・rump:
演算精度の桁数（デフォルト1）:

ターミナルで

$ pycalcl

で過去入力コードを記録した
$HOME/myscript/pycalcl.txt
を開く

【5通りの出力】

pprint(eq,use_unicode=False)
pprint(eq,use_unicode=True)
pprint(eval('eq'))
print(eval('eq'))
pprint(N(eq,{N:- 1}))

【インストール】

HOMEにディレクトリmyscriptをつくる
HOME/myscriptに
pycalc.sh
pycalcl.txt
を配置

MacOSのzshの場合.zshrcに
ubuntuのBashの場合.bashrcに
に以下を追記

source $HOME/myscript/pycalc.sh

一行数式SymPy計算機シェルスクリプト pycalc.sh

# pycalc.sh
function pycalcl(){
open HOME/myscript/pycalcl.txt
}

function pycalc(){
echo '【One liner SymPy Calculator（一行数式SymPy計算機）】Ver.20230508'
echo '　商// 剰余% 商と剰余divmod()分数2/3 Frac(2,3) 平方根sqrt(2) 三角関数sin(pi/3) 指数exp() 自然対数log(E**2) 常用対数log(2,10)'
echo '　虚数j(標準) (2+3j)*(5-7j)'
echo '　虚数I(SymPy) exp(cos(E**I)+sin(E*pi)) I**I'
echo '　素因数分解 factorint(1000) 階乗factorial(10)'
echo '　代数演算 シンボリック変数aからz expand((x+y)**10) factor(a**10-b**10) '
echo '　数列 Sum(k**2,(k,1,n)).doit()'
echo '　方程式 solve(a*x**2+b*x+c,x) 連立方程式 solve([x+y-4,x-y-2],[x,y])'
echo '　微分方程式 変数はf限定 dsolve(Eq(f(t).diff(t, t) - f(t), exp(t)), f(t))'
echo '　微分 diff(x**2,x) 積分integrate(x**3,x) 定積分integrate(x**3,(x,0,1))'
echo '　無限oo integrate(1/(1+x**2), (x, -oo, oo)) '
echo '　テイラー展開 series(sin(x),x, 0, 12)'
echo '　テイラー展開 係数リスト taylor(sin, 0, 5)'
echo '　行列 Matrix([[1, 2], [2, 2]]).eigenvals()'
echo '　関・ベルヌーイ数 bernoulli()'
echo '　ゼータ zeta() zetazero()'
echo '　ブール値検算 1+1 == 3   expand((x+y)**2) == x**2 + 2*x*y + y**2'
echo '　グラフ plot(x**2, (x, -1, 2), ylabel = "y")'
echo '　mpmath任意精度浮動小数点演算パッケージによる精度計算 (-2)**mpf("0.5")'
echo '　Rumpの例題 入力にrump https://ja.wikipedia.org/wiki/%E7%B2%BE%E5%BA%A6%E4%BF%9D%E8%A8%BC%E4%BB%98%E3%81%8D%E6%95%B0%E5%80%A4%E8%A8%88%E7%AE%97'

echo -e '\n数値計算式・代数式（変数eq）・rump:'
read eq
echo -n '\n演算精度の桁数（デフォルト1）:'
read N
echo "eq" >> HOME/myscript/pycalcl.txt
echo -E "
from mpmath import *
mp.pretty = True
from sympy import *
from fractions import Fraction as Frac
init_printing()
var('a:z')
f = Function('f')

rump = 'rump'
ifeq == 'rump':
    a=77617
    b=33096
    c=333.75*b**6+a**2*(11*a**2*b**2-b**6-121*b**4-2)+5.5*b**8+a/(2*b)
    print(f'通常計算{c}')

    def g(a, b):
        return (mpf('333.75')*b**6 + a**2*(11*a**2*b**2-b**6-121*b**4-2)+mpf('5.5')*b**8+a/(mpf('2')*b))
    for mp.dps in range(1, N+1):
        print(mp.dps,g(mpf('77617'), mpf('33096')))
else:
    mp.dps ={N:- 1}
    # mp.prec = 200
    print('\npprint(eq,use_unicode=False)')
    pprint(eq,use_unicode=False)

    print('\npprint(eq,use_unicode=True)')
    pprint(eq,use_unicode=True)

    print('\npprint(eval(eq))')
    pprint(eval('eq'))

    print('\nprint(eval(eq))')
    print(eval('eq'))

    print('\nN(eq,{N:- 1})')
    pprint(N(eq,{N:- 1}))
" | python
}

使用例

初等関数・分数

スクリーンショット 2023 05 08 12 00 14

虚数

スクリーンショット 2023 05 08 12 00 35

分数

スクリーンショット 2023 05 08 12 00 50

方程式

スクリーンショット 2023 05 08 12 01 05

逆行列

スクリーンショット 2023 05 08 12 01 19

微分方程式

スクリーンショット 2023 05 08 12 01 32

検算

スクリーンショット 2023 05 08 12 01 46

検算

スクリーンショット 2023 05 08 12 01 59

任意精度浮動小数点演算

スクリーンショット 2023 05 08 12 02 43

Rumpの例題の検証

スクリーンショット 2023 05 08 12 03 21

投稿日:

関・ベルヌーイ数　Pythonコード

■12月12日(日)11:00-12:00 桜井進のPython･UNIX・Math教室（応用コース）
スクリーンショット 2021 12 13 1 20 37

超入門・Pythonで関・ベルヌーイ数

  1. My code
  2. Akiyama-Tanigawa algorithm
  3. B(n) is just sum of k^n formula linear term coefficient.
  4. Zeta function algorithm

Seki bernoulli 001
Seki bernoulli 002
Seki bernoulli 003

# 1. My code
# seki-bernoulli_1.py

from fractions import Fraction

# 二項係数 Combination nCk
def comb(n, k):
    prod = 1
    for t in range(min(k, n-k)):
        prod = prod * (n-t)//(t+1)
    return prod

# Seki-Bernoulli number Bn
def B(n):
    if n == 0:
        return 1
    else:
        ss = 0
        for k in range(0, n):
            ss = ss + ((-1)**k) * comb(n+1, k)*B(k)
        return Fraction((-1)**(n+1), n+1) * ss


n = int(input('B(n)のnの上限 >>> '))
import sympy
for i in range(n + 1):
    print(f"B({i})= {B(i)}".ljust(20, " "),f"sympy.bernoulli({i})=",sympy.bernoulli(i))

実行結果

$ python seki-bernoulli_1.py
B(n)のnの上限 >>> 20
B(0)= 1              sympy.bernoulli(0)= 1
B(1)= 1/2            sympy.bernoulli(1)= -1/2
B(2)= 1/6            sympy.bernoulli(2)= 1/6
B(3)= 0              sympy.bernoulli(3)= 0
B(4)= -1/30          sympy.bernoulli(4)= -1/30
B(5)= 0              sympy.bernoulli(5)= 0
B(6)= 1/42           sympy.bernoulli(6)= 1/42
B(7)= 0              sympy.bernoulli(7)= 0
B(8)= -1/30          sympy.bernoulli(8)= -1/30
B(9)= 0              sympy.bernoulli(9)= 0
B(10)= 5/66          sympy.bernoulli(10)= 5/66
B(11)= 0             sympy.bernoulli(11)= 0
B(12)= -691/2730     sympy.bernoulli(12)= -691/2730
B(13)= 0             sympy.bernoulli(13)= 0
B(14)= 7/6           sympy.bernoulli(14)= 7/6
B(15)= 0             sympy.bernoulli(15)= 0
B(16)= -3617/510     sympy.bernoulli(16)= -3617/510
B(17)= 0             sympy.bernoulli(17)= 0
B(18)= 43867/798     sympy.bernoulli(18)= 43867/798
B(19)= 0             sympy.bernoulli(19)= 0
B(20)= -174611/330   sympy.bernoulli(20)= -174611/330

案の定、関数B(n)は再帰的定義をしているため、B(30)の計算はできません。
結果をsympy.bernoulli()と比較してみます。

Seki bernoulli 002

このアルゴリズムはヤコブ・ベルヌーイによるべき乗和による導入によるもの
したがって、B(1)=\frac{1}{2}
sympy.bernoulli(1)=-1/2 は、現在主流である母関数\frac{x}{e^x-1}による関・ベルヌーイ数の定義

Seki bernoulli 001

次はゼータの負の整数に対する公式によるもの
関数mpmath.zeta()を用いるだけのコード

# 4 zeta fuction algorithm
# seki-bernoulli_4.py

print('4 zeta fuction algorithm')
n = int(input('B(n)のnの上限 >>> '))

from mpmath import zeta
def B(m):
    b = -m * zeta(1-m)
    return b

from sympy import *
for i in range(1, n + 1):
    b = B(i)
    c = Rational(b)
    print(f"B({i})= {c}".ljust(20, " "),f"sympy.bernoulli({i})=",bernoulli(i))

実行結果

$ python seki-bernoulli_4.py
4 zeta fuction algorithm
B(n)のnの上限 >>> 50
B(1)= 1/2            sympy.bernoulli(1)= -1/2
B(2)= 1/6            sympy.bernoulli(2)= 1/6
B(3)= 0              sympy.bernoulli(3)= 0
B(4)= -1/30          sympy.bernoulli(4)= -1/30
B(5)= 0              sympy.bernoulli(5)= 0
B(6)= 1/42           sympy.bernoulli(6)= 1/42
B(7)= 0              sympy.bernoulli(7)= 0
B(8)= -1/30          sympy.bernoulli(8)= -1/30
B(9)= 0              sympy.bernoulli(9)= 0
B(10)= 5/66          sympy.bernoulli(10)= 5/66
B(11)= 0             sympy.bernoulli(11)= 0
B(12)= -691/2730     sympy.bernoulli(12)= -691/2730
B(13)= 0             sympy.bernoulli(13)= 0
B(14)= 7/6           sympy.bernoulli(14)= 7/6
B(15)= 0             sympy.bernoulli(15)= 0
B(16)= -3617/510     sympy.bernoulli(16)= -3617/510
B(17)= 0             sympy.bernoulli(17)= 0
B(18)= 43867/798     sympy.bernoulli(18)= 43867/798
B(19)= 0             sympy.bernoulli(19)= 0
B(20)= -174611/330   sympy.bernoulli(20)= -174611/330
B(21)= 0             sympy.bernoulli(21)= 0
B(22)= 854513/138    sympy.bernoulli(22)= 854513/138
B(23)= 0             sympy.bernoulli(23)= 0
B(24)= -236364091/2730 sympy.bernoulli(24)= -236364091/2730
B(25)= 0             sympy.bernoulli(25)= 0
B(26)= 8553103/6     sympy.bernoulli(26)= 8553103/6
B(27)= 0             sympy.bernoulli(27)= 0
B(28)= -23749461029/870 sympy.bernoulli(28)= -23749461029/870
B(29)= 0             sympy.bernoulli(29)= 0
B(30)= 8615841276005/14322 sympy.bernoulli(30)= 8615841276005/14322
B(31)= 0             sympy.bernoulli(31)= 0
B(32)= -7709321041217/510 sympy.bernoulli(32)= -7709321041217/510
B(33)= 0             sympy.bernoulli(33)= 0
B(34)= 2577687858367/6 sympy.bernoulli(34)= 2577687858367/6
B(35)= 0             sympy.bernoulli(35)= 0
B(36)= -26315271553053477373/1919190 sympy.bernoulli(36)= -26315271553053477373/1919190
B(37)= 0             sympy.bernoulli(37)= 0
B(38)= 2929993913841559/6 sympy.bernoulli(38)= 2929993913841559/6
B(39)= 0             sympy.bernoulli(39)= 0
B(40)= -261082718496449122051/13530 sympy.bernoulli(40)= -261082718496449122051/13530
B(41)= 0             sympy.bernoulli(41)= 0
B(42)= 1520097643918070802691/1806 sympy.bernoulli(42)= 1520097643918070802691/1806
B(43)= 0             sympy.bernoulli(43)= 0
B(44)= -27833269579301024235023/690 sympy.bernoulli(44)= -27833269579301024235023/690
B(45)= 0             sympy.bernoulli(45)= 0
B(46)= 596451111593912163277961/282 sympy.bernoulli(46)= 596451111593912163277961/282
B(47)= 0             sympy.bernoulli(47)= 0
B(48)= -5609403368997817686249127547/46410 sympy.bernoulli(48)= -5609403368997817686249127547/46410
B(49)= 0             sympy.bernoulli(49)= 0
B(50)= 495057205241079648212477525/66 sympy.bernoulli(50)= 495057205241079648212477525/66
投稿日:

解決問題「コラッツ・角谷予想」（後半）

空間情報連載
Pythonで数学を学ぼう！

桜井進のPython＋Math教室 第9回

Pythonで数論！未解決問題「コラッツ・角谷予想」（後半）
1.コラッツ予想

2.コラッツ・シークエンスのステップ数
コラッツ・シークエンスのステップ数を求めるコード

# colg.py
# ステップ数（最大ステップ数）と頻度のグラフ描画
import numpy as np                  # 配列を扱う数値計算ライブラリNumPy
import matplotlib.pyplot as plt     # グラフ描画ライブラリmatplotlib
import japanize_matplotlib          # matplotlibの日本語化
import datetime as dt
from decimal import Decimal

while(1):
  Model = input('1.a≦n≦bに対するコラッツ・シークエンスのステップ数の最大値とnを算出\r\n'
                '2.a≦n≦bに対する横軸ステップ数、縦軸頻度の棒グラフ描画\r\n'
                '1、2のどれかを入力 ')
  if Model.isdecimal():
    Model = int(Model)
    if 1 <= Model <= 2:
      break