お年賀 2022年

さて、久しぶりにHP Primeのシミュレータで動くプログラムを作業してみました。
少々どころか、かなり長いコードなので、HP Prime実機のキーをポツポツ叩いて打ち込むのは苦行です。
PCとかを経由して、手軽に実機へコードを送り込む方法とかあるのかな ?
当方、未だに実機を入手しておらんので、この辺りの事情を御存知の方は、お知らせ戴けると幸甚であります。

実をいいますと、HP Primeでpythonが動く、と書いた手前、pythonで動くプログラムをでっち上げようと思ったのですが、global 文が機能しない様な感触だった (古いシミュレータだから ?)ので、今回はPPLで書き下した次第。

追記 on 2022-03-11

別の所に、画像を変更したPPLコードを掲載しました。
https://hpcalc-fun.blogspot.com/2022/03/2022.html

// title : graphic example
// begin : 2022-01-04 04:10:54 

//
EXPORT colour;

//  graphics print func
gprint(I, J, str) 
BEGIN
LOCAL c;
  D := 4;

  FOR X FROM 1 TO DIM(str) DO
      c := MID(str, X, 1) ;
      N := INSTRING("0123456789ABCDEF", c) ;
      L := colour(N) ;
      RECT_P(I, J, I+D, J+D, L) ;
      I := I+D ;
  END;
END;

//  
EXPORT main()
BEGIN

RECT();

//  palette
colour := {} ;
colour := CONCAT(colour, RGB (32, 30, 22) );
colour := CONCAT(colour, RGB (93, 94, 86) );
colour := CONCAT(colour, RGB (51, 69, 71) );
colour := CONCAT(colour, RGB (77, 52, 30) );
colour := CONCAT(colour, RGB (111, 146, 156) );
colour := CONCAT(colour, RGB (94, 121, 131) );
colour := CONCAT(colour, RGB (173, 199, 205) );
colour := CONCAT(colour, RGB (228, 229, 222) );
colour := CONCAT(colour, RGB (163, 160, 155) );
colour := CONCAT(colour, RGB (172, 138, 110) );
colour := CONCAT(colour, RGB (153, 113, 87) );
colour := CONCAT(colour, RGB (210, 173, 158) );
colour := CONCAT(colour, RGB (200, 150, 119) );
colour := CONCAT(colour, RGB (223, 201, 180) );
colour := CONCAT(colour, RGB (98, 76, 50) );
colour := CONCAT(colour, RGB (144, 93, 56) );

// ======== graphic print

gprint(  0 ,  0 ,  "000000012324445205444467526677846777425466444541122200214444467776645000445"  ) ;
gprint(  0 ,  4 ,  "000000084114445505644467226777846666405666444555122000214114667776445002465"  ) ;
gprint(  0 ,  8 ,  "00000091A845444504444466524777446667405467454555122000214444777677445002445"  ) ;
gprint(  0 ,  12 ,  "000000432664466525645466554777446676404476456444112202111118667776442002445"  ) ;
gprint(  0 ,  16 ,  "00000012046676640544446752667784666760546664688B88BC88DB6954777776442002645"  ) ;
gprint(  0 ,  20 ,  "00100012246677640564446742477744666762546768D6DB8BBBBCBDDDB8667777642002645"  ) ;
gprint(  0 ,  24 ,  "0022011221666774254444664247776566666546D698DBD77DDD77DD8D78886777645002445"  ) ;
gprint(  0 ,  28 ,  "00200522224677642444446742477764667767777BEE8BB77D777777DBB7644666445002465"  ) ;
gprint(  0 ,  32 ,  "2332252322567764054444674267776467777777745399D777D777777DB8D68488641202645"  ) ;
gprint(  0 ,  36 ,  "3E2025002E5467742244446742467784677777771B828BD7D77D7777DDD8968886684002442"  ) ;
gprint(  0 ,  40 ,  "2E30220302114774226444676267776BBD77777B8123BDB67D7777777DD8589487766502645"  ) ;
gprint(  0 ,  44 ,  "322312020222477425444466424777BBBDD7777B8308BBD7DDDD7777DDC9911EAD777B55445"  ) ;
gprint(  0 ,  48 ,  "323212030112567422444466424777BBD77777D89EA8BDDBDDDDD7D7DB9891AA99B77768645"  ) ;
gprint(  0 ,  52 ,  "2302100222242674225555468256DBBB77777777B86BBBBDDD6BDDDDB8888596D899D776645"  ) ;
gprint(  0 ,  56 ,  "3222100202565464224444666058B7BBBD777777DBDBDBBDBBBBBBBB8888888D77D89877645"  ) ;
gprint(  0 ,  60 ,  "E32110020E264266224554676286BBBBDBD77BDD8BBD6BD6DBB88BB88CAA888D777D8967745"  ) ;
gprint(  0 ,  64 ,  "E222E0020E16826622445546B4BA8BBBBB9BBCB98BBBBBBBBB8B8B989BA195BD7777789D7D6"  ) ;
gprint(  0 ,  68 ,  "3222200E0228656620445466BBAABBBBBBA9AA89BBBB6BBBD88BBBB989AA9AB7DD777DB8777"  ) ;
gprint(  0 ,  72 ,  "202010020114744422445866B9ACBBBBBCAAB8B89BBBBBA8BBBBB8BCC9AA9ABDDDD7777B877"  ) ;
gprint(  0 ,  76 ,  "E22210020218764422444888BCACBBBBBCABBBB8B595BB9188BBBB888A19A8BB98DD7777B8D"  ) ;
gprint(  0 ,  80 ,  "E20212020E28765422444DBBBBAABBBDBBBBDBBBCA851BBAA8B8BB8C9AA9A8BAA9DDD7777DD"  ) ;
gprint(  0 ,  84 ,  "2022E0000228774550448DBBBBAABBBBBBBBBBBB91B82ABB1BB8BC9C9A99199AF989D77777D"  ) ;
gprint(  0 ,  88 ,  "22321002211477655245BDDBBBCCABBBBBBBBBBBB14129B8F98B8B989A9F9AAAE999DDDDD77"  ) ;
gprint(  0 ,  92 ,  "222220020114676450448DDBBBCCCABDBBBBBBBB812238BBAA8B99CA999E91A1FA99B9DD677"  ) ;
gprint(  0 ,  96 ,  "222EE00301247774525486DBDBBCCABBDDDBDBBBB5103BBBAAB8999CA1AAEEA11AA9C9CD8DD"  ) ;
gprint(  0 ,  100 ,  "323222002229777652445DDBBBBBCBBBDBBDBDBBBA1E1B8BAA9CC999A9AAEEAE1AA9A988DDD"  ) ;
gprint(  0 ,  104 ,  "212220020E156776554448BBBBDDBDDDDDBBBBBBB95AABBB998C999119A1EF1E1FA9F9988DB"  ) ;
gprint(  0 ,  108 ,  "1112200522E577774446488DDDDBBDBBBBDBBBBBBCB9BCCC9C9999AEF91EE11E11AAE999CC9"  ) ;
gprint(  0 ,  112 ,  "21125200022167674445466BBDDBDDDBBBBBBBBBCCCCBBB89CC99AEE19E0E1E1EFAE1A998B8"  ) ;
gprint(  0 ,  116 ,  "151E12000144666654812246BDBDDBDBBBBBBBBCBCC8CCC998CAAFEEAAE3EEE11FE21A599CC"  ) ;
gprint(  0 ,  120 ,  "151222000224666686DD7641CDBDDDBDBBBCCBCCB99CCB8C8CAAEE3191321EE1FEEEFAA998C"  ) ;
gprint(  0 ,  124 ,  "1122E200000224686DD77777ABBBDBBBBBBBCCCC99C89C9CCAFEEE319E33E1FE1E2E11AA9C8"  ) ;
gprint(  0 ,  128 ,  "21E2E2000112039DDDD7777DCCCDBBDBBBC8CC9C99C9C8CC9F133E3AA332EEE1EEEEEEAAA99"  ) ;
gprint(  0 ,  132 ,  "212222200111199CDD77D77DCCBBBDDBBBBCC9C99C9C9CBAA133E3EA323EEEEEEEE1EF1A199"  ) ;
gprint(  0 ,  136 ,  "E22E2E200221999B7D7D77DCBCCDBBBBBB9999AA9AA999A1E33333AF0333EEE1EEE1EFAA199"  ) ;
gprint(  0 ,  140 ,  "11223200022289CD77DD77CCCCCDBBBCCCC9999C9999A12E33003FA003E22E1EEE1E211AF99"  ) ;
gprint(  0 ,  144 ,  "222322200221CCDDDDDDDDDCCCCBBBDBB99A98C9A9AEE3030300FA0033E23E1EEE1E2F1F199"  ) ;
gprint(  0 ,  148 ,  "32020230022ACDDDCCBBCCCACCCCCDBCCCAAAA9AAA1EF3300303F30330333113EEE3E11A19A"  ) ;
gprint(  0 ,  152 ,  "022002000EABDDD9FCCACFFFCCBCBCCCCAAAF1A58AAAA133303F300330322EE3EEE22EA1F9A"  ) ;
gprint(  0 ,  156 ,  "302000001B7DDDCCCCAC9FFFFCCCCCCAAAAFAF1ABAB33A1103E300303203E1332EE2E1AE19A"  ) ;
gprint(  0 ,  160 ,  "00000018DDDDDDC9CFFFFFEEF33EF3FFFFEFEEEEA3193FB9F330033003333303EEE23EF1FA1"  ) ;
gprint(  0 ,  164 ,  "00001BDDDDDC9CA9FFFFEEE3333333EE1FA1FFAFA1F13EFA3333333303030003E1E0EEA119A"  ) ;
gprint(  0 ,  168 ,  "0E9DDD7DBCFCA9FFFFEF333FE33FF9AFAFE33333EF1E1E33E333303003030003EA10EE1EA9A"  ) ;
gprint(  0 ,  172 ,  "8DD7DDCCFCFAFFFFFEE33FF33FAA9FFF3303330300003F11EE33300000000003FAE03EEAF9A"  ) ;
gprint(  0 ,  176 ,  "DDCCCFCFACFFFEEFFEFFFFFFCCFFFFE303FCFFE3000000003000000000000000AAE2EE1F19F"  ) ;
gprint(  0 ,  180 ,  "CFCCAACAAFFFFFFF3FFA9CCFA33F9F33F9CAAFF3300000000000000000000000FA122EF1F9A"  ) ;
gprint(  0 ,  184 ,  "FCFFCFFFFEFFFFFFACAAFFFAF3FCA33ACCAAAFFF330000000000000000000000E99E2EF119A"  ) ;
gprint(  0 ,  188 ,  "FFCFFFFFEFFFAFCCFFFEFCCA3FCCF3ACCAACAAFFF33000000000000000000000E98EEE11AAA"  ) ;
gprint(  0 ,  192 ,  "CCFFFEF3FFFFCCFFFFACCCAFFCCF3ACC9CCAAC9BAFF300000000000000000000EFCAF3A3F9A"  ) ;
gprint(  0 ,  196 ,  "FFFFFFFFFFCCFFFFACCCCCF3FC9FFCBCCCCBBDD77DBA30000000000000000000EFC9EEFEA99"  ) ;

//  show screen and wait
FREEZE;

END;

python電卓花盛り

永い事、「買いたいなぁ」と思っている内に、HP Prime は軒並み値上がり、こうなると、金欠の当方としては、なかなか手が出ない。
追い打ちを掛ける様に、いよいよ、TI-84 Plus CE pythonが「アメリカからの特別のプログラムを送ります」という感じで、出回り始めました。
目下、python電卓が豊作となっており、どれを買ったらいいもんじゃろか ?

そんな年の瀬、CASIO fx-CG50のアプデが提供されておりますネ。

https://support.casio.jp/download.php?cid=004&pid=1067

「オッ、コレは !?」と思っておりました。藤堂様、来年と言わず「Buy Now !」なのかも ?

https://twitter.com/ShunsukeTodo/status/1463412223888420865

しかし、今回のアプデは、pythonの拡張ではなく、統計分布アプリの新設、というものでした。
まあ、今のpython機能でも十分なものはあると思いますが、徐々にでも拡張がされたらウレシイかなぁ。

五輪の夏、DELTAの夏 - 夏休みの自由研究 2021年

目下、新型感染症の猛威は留まる所を知りませんが、今日流行しているものは「デルタ株」という新しい型らしい。

ref. 埼玉県 インドで確認「L452R」変異ウイルス 陽性率85.4％に - NHK news https://www3.nhk.or.jp/news/html/20210811/k10013195471000.html

東京都は、新型肺炎の感染についての情報を公開しており、 中には、デルタ株についての検査報告もあがっております。

ref. https://github.com/tokyo-metropolitan-gov/covid19/blob/development/data/variants.json

 

日付	デルタ株件数	全件数
2021-04-30	0	76
2021-05-03	1	121
2021-05-10	2	103
2021-05-17	8	139
2021-05-24	3	372
2021-05-31	15	309
2021-06-07	32	1002
2021-06-14	127	1516
2021-06-21	261	1770
2021-06-28	502	2336
2021-07-05	934	3050
2021-07-12	1948	4220
2021-07-19	3674	5688
2021-07-26	9983	12222

このデータを使い「不要不急の外出」を控えつつ「夏休みの自由研究」をしてみましょう。

上記の情報ですが、検査数が一定していないので、全体の感染数に対するデルタ株の割合を調べてみます。 グラフ電卓の散布図プロット機能を使うと、デルタ株の割合の時系列変化がよく判ります。 デルタ株の割合のグラフを描くと、どうも「ロジスティック関数」のようであります。

ref. ロジスティック方程式 - Wikipedia https://ja.wikipedia.org/wiki/%E3%83%AD%E3%82%B8%E3%82%B9%E3%83%86%E3%82%A3%E3%83%83%E3%82%AF%E6%96%B9%E7%A8%8B%E5%BC%8F

最近のグラフ電卓には「ロジスティック回帰」の機能がありますから、これで調べてみるとよさそうですが、古いものには実装されておりません。

そこで、古い電卓でも分析できる様に、作業してみました。

ロジスティック方程式ですが、「全体の割合」という事で限定すると、つぎの式で記述できます。

これを書き換えると、

となって、指数回帰の機能で回帰分析出来そうです。

実際にロジスティック回帰を実施するには、日付を週番として独立変数x に、「デルタ株」件数/全件数を従属変数yとしますが、従属変数yを(1-y)/yと変数変換して回帰分析を実行します。

デルタ株件数が0になっている所では、この変数変換が計算できませんから、ごく小さな値を設置して、計算エラーを回避するという工夫が必要になります。

こうして指数回帰を実施すると、つぎの式が得られました。

2101.37490861*0.480788843194^X 

もちろん、この式は指数回帰の式でありますから、これをロジスティック回帰の式に書き換えると、

1/(1+2101.37490861*0.480788843194^X)

が得られます。

グラフを描くと、こんな具合になります。

ここでは、HP Prime シミュレータを利用しましたが、HP Prime には、ロジスティック回帰の機能がありますネ。 しかし、ロジスティック回帰に特化した計算方式ではないらしく、上記のデータで実行しても、うまいこと計算してくれません。ｼｸｼｸ。 そこで、今回のような姑息な手段を講じたのですが、お手持ちの高機能電卓でお試し戴く際の参考にして戴きたく思います。

デルタ株、猛威をふるっており、全国的にも感染者数が激増しております。

接種をしたからといっても、まだまだマスクは手放せない模様。くれぐれもお気をつけて頂きます様。

 

HP Prime CASの計算の数値範囲

HP PrimeのCASは色々と楽しく、ネットで見掛けた問題を幾つか解くのに使ってみました。

1. x^sqrt(x) = x*sqrt(x)

Can you solve this rational equation ?

cf. https://www.t-3.com/thinking/brilliant-orgs-brilliant-crm/

CASで

fsolve(x^sqrt(x)=x*sqrt(x), x)

と入力します。すると、一度、こんな警告が出ますが、

ENTERキーを押して送ると、[1, 2,25] と、答えを返します。ウマイッ !

2. x! = (7!)! / 7!

what is x ?

cf. https://sites.google.com/a/g.coppellisd.com/coppell-ib-math/math-sl/topics---sl/probability/counting-principles

これは、すんなりと解けるものではなく、別に色々とやって答が判明していたので、その検証に使おうとしました。

答は x=5039 となるらしい。PCで動くmaximaで手当たりしだいに試し結果を得たですが、検証をHP PrimeのCASで出来ないか、試した次第。

(7!)!/7!-5039!
    undef

所が ... 計算できまへん。undefとなってしまいます、アカン。

では、何処までの階乗が計算できるのか、調べてみました。どうも、1006! までは計算できますが、1007! になるとお手上げです。ウーム。

ちなみに、1006! がどんな数値なのか、そこを見ておこうと思います。
整数では結果が得られますが、結構桁数があります。浮動小数点表記だと、どうなるのか ?

HP Primeの浮動小数点表示は、10 E 500 以上でエラー、もしくは Inf となります。1006! は、この範囲にはないでしょう。そこで、対数を取って階乗を計算し、それを浮動小数点表記にする、という作業を行います。

log(1006!) の計算

ΣLIST(MAKELIST(LOG(X),X,1,1006,1))  Shift+ENTER (概数計算の入力)
    2585.61374471

これより、1006! = 4.10908107239 E 2585 と概算されます。
(10^0.61374471 = 4.10908107239 ですネ)

log(1007!) の計算

ΣLIST(MAKELIST(LOG(X),X,1,1007,1))  Shift+ENTER
    2588.616774418

これより、1007! = 4.13784463462 E 2588 と概算されます。

浮動小数点数は10E500まででしたが、整数は、この程度まで計算してくれる模様。
そこで、限界を調べるべく、電卓のキー (実際にはシミュレータのボタン)を叩きまくり、つぎの様な値までは扱えるらしい事が判りました。

3605227827*10^2579 = 3.605227827 E 2588

そして、この整数は、おそらく2進数integerで内部的に保持しているであろうから、何ビットで表現できるのか、と考えてみました。

2を底とする対数を取ったりしてみると、2^8599 が3.605 E 2588 に近い値となる様子。
HP Prime CASでは、2^8599以上の整数値は扱えないらしい、と考えられます。

所が、こんな結果が !

s=2^8598
    1802613913 ... 2262841344 (途中省略の表記です)
s*s
    undef

ここまでは当然ですネ。

s=2^8597
    9013069568 ... 1131420672
s+s
    1802613913 ... 2262841344
s*4
    undef

2進数的に1桁落として加算します。ここも十分理解できます。

s+s+s+s
    3605227828 ... 4525682688

ムムッ !? コレはオカシイ、undefにならんじゃないの !!

調子に乗って、64回加算を行いましたが、一応結果が得られました。
加算の処理が簡易だからなのか、回数を増やしても、ある程度ならば計算してくれる様です。

こんな具合で、整数の場合、3.6 E 2588 程度まで計算してくれる事が判りました。
欲を書くと、もっと長い桁数の整数計算が出来ると良かったのですが、普通の関数電卓では、ここまで扱えないのが当たり前なので、コレ以上の計算の需要はPCの利用を、という事ですネ、ハイ。

HP Prime CAS でやってみた

こんな話題をみつけました。

cf. https://twitter.com/dentakuin/status/1145548793527459840

陰関数表示のグラフ描図との事です。

以下のコードは、sin(x^2)+sin(y^2)>=1 を描く、CASのpythonコードの例です。


#cas

def sinsqs() :
  wi = 320
  he = 240
 
  Xmin = -6.4
  Xmax =  6.4
  Ymin = -4.8
  Ymax =  4.8
 
  for k in range(he) :
    for j in range(wi) :
      x = Xmin + (Xmax-Xmin)*j/wi
      y = Ymax - (Ymax-Ymin)*k/he
      if sin(x**2)+sin(y**2) >= 1 :
        col = RGB(255,0,0)
      else :
        col = RGB(255,255,255)
      PIXON_P(j,k,col)
  WAIT

#end

HP PrimeのCASで円周率の計算

HP PrimeのCASで円周率の計算をやってみました。

pythonでは多桁整数の計算が出来ますが、PrimeのCAS・XCASでも同様なので、これを利用してMachinの公式による円周率計算を行った次第。

つぎのプログラムをHP Primeのプログラムに書き込み、実行後、CASに移動して

m(100)

の様に打ち込んでやります。すると、円周率の計算を行って、結果が出てきます。

項の計算で末尾まで行っているわけではないので、100と入力しても、100桁の精度があるのかはわかりません。100桁の精度が所望の場合には、120とか、大きな値を設定して実行するとよいと思います。
また、このプログラムもモノホンのpythonで実行できると思います。

#cas
def  machin(n) :

  s = 0
  t1 = 10**n // 5
  t2 = 10**n // 239

  f = 0

  for i in range(n):
    d = 1+2*i
    if f == 0 :
      s = s + t1//d * 4
      s = s - t2//d
      t1 = t1 // 5 // 5
      t2 = t2 // 239 // 239
    else :
      s = s - t1//d * 4
      s = s + t2//d
      t1 = t1 // 5 // 5
      t2 = t2 // 239 // 239
    f = 1 - f

  print(s*4)

#end

HP PrimeのCASでお遊び

HP PrimeのCASは「XCAS」を使っているそうで、独自スクリプトのほかに、python文法に準じたスクリプトを実行できる、との事。
そこで、簡単なスクリプトを作成してみましたヨ。

参考にしたのは、つぎのTweet。

cf. https://twitter.com/RR_Inyo/status/1217489174128758784

少々古いものですが、色々と研究して、どうにか再現できた次第。

使い方

1. HP Primeのプログラム・エディタで、適当な名前でファイルを作成します。この時、CASのチェックボックスにチェックを入れるのを忘れずに。
2. プログラムを実行します。しかし、この時点では、何も実行されません。
3. CASを呼び出し、iter_cos(20) と入力して実行します。

コード解説

HP PrimeのXCASではPPLのグラフィクス命令が使えるので、こんなコードが動くのですネ。
また、XCASのコードなので、頭と尻に#CAS, #END というタグが付いておりますが、同時に、ファイル名として準備されるスクリプトヘッダは、必要無さそうな感じなので削除しております。
スクリプトを自動実行できない所はありますが、グラフィクスのコードも動くので、そこそこ楽しめそうです。

コード

#cas

def iter_cos(n):
  width = 320
  height = 240
  for k in range(height):
    for j in range(width):
      x = j/width*10.0-5.0
      y = 5.0-k/height*10.0
      z = x+y*i
      l = 0
      while l<n and abs(z)<50:
        z = COS(z)
        l=l+1
      col = rgb(l/20*255,0,0)
      PIXON_P(j, k, col)
  WAIT

#end

2週間遅れの「バレンタインデー」

先日、こんな情報をみたのでした。

https://twitter.com/NumWorks/status/1360951932970409987

話によると、つぎの様な数式でグラフが描けるのだとか。

  (sqrt(cos(x))*cos(400*x)+sqrt (abs(x))-0.4)*(4-x*x)^0.1

cf. Google ‘Heart Graph’: Geeky Surprise Will Bring You Joy On Valentine’s Day (PICTURES, VIDEO)
https://www.huffpost.com/entry/google-heart-graph-valentines-day_n_1276391

で、HP Primeシミュレータでやってみましたよ、久しぶりに !

そこそこの絵が出て来ます。一部、塗りつぶしがうまく行かないのですが、ズームで細部を確認戴くと、理由が判るものと思います。

数式の (4-x*x)^0.1 の部分から、値域 {x | -2 < x < 2} で実数値になりますネ。
また、cos(400*x) の部分で塗りつぶしを実現している様です。

しかし、よくもまぁ、こんな数式があるもんです。驚いちゃったよ。

今日はホワイトデーなので (?)、久しぶりにHP Primeのコード

まずは、こちらのtweetを御参照

cf,
https://twitter.com/RR_Inyo/status/1217489174128758784

z' = cos(Z) の繰り返しで、こんなグラフィクスが描けるそうで、HP PRIMEでもやってみましたヨ。
まだ、実機を持っていないので、今回もSIMで作業ですが、以下のコードを御利用下さい。

// title : iter_cos.ppl - iterative complex cosine graph for HP Prime
// begin : 2020-02-16 20:37:33 

//
pixVal(z, c)
BEGIN
  LOCAL  cnt := 0;

  //  この判断式がミソです
  //  Mandelobrot, julia setでは、abs(z)<4 .0="" br="">  //  このグラフィクスのバヤイはこんな風でした
  WHILE (ABS(z) <= 50) AND (cnt < 20) DO
    // modify this form to enjoy
    z := COS(z);
    cnt := cnt+1;
  END;

  //  color value : 0 for black, 16777215 for white
  RETURN RGB(1*cnt/20*255, 0.75*cnt/20*255, 0.25*cnt/20*255);

END;

EXPORT itercos()
BEGIN
  // Clean the screen (G0):RECT();

  LOCAL dx, dy, z, xp, yp; 
  LOCAL xmin, xmax, ymin, ymax;
  LOCAL pixStat;

  //  screen range
  xmin := -10;
  xmax :=  10;
  ymin :=  -5;
  ymax :=   5;

  dx := (xmax-xmin)/320;
  dy := (ymax-ymin)/240;

  // we loop over each pixel
  // Of the HP Ptime screen 
  FOR yp FROM 0 TO 239 DO
    FOR xp FROM 0 TO 319 DO
      z := (xmin+xp*dx, ymax-yp*dy);
      pixStat := pixVal(z, Z0);
      PIXON_P(xp, yp, pixStat);

    END;
  END;

  //  key wait loop
  REPEAT UNTIL GETKEY() == -1;
  FREEZE;

END;

Numworks, mandelbrotのコードを改変したとの事で、単純に式を置き換えたのですが、それだけでは、こうしたグラフィクスにはなりませんでした。
色々といじくり回し、最終的に繰り返しの判断部分をいじる事で、Tweetの様なグラフィクスを得た次第。

最近、HP Prime G2の価格が大分下落してきております。

cf.
https://www.amazon.co.jp/gp/offer-listing/B07HF6RXGG/ref=dp_olp_new?ie=UTF8&condition=new

G1との価格差も狭まってきました。Sentaro様曰く、G2は動作軽快でメモリ増量。ウーム。

一方、G2は、新しいファームがあるそうです。

cf. HP Prime Beta Software
https://www.hpcalc.org/prime/beta/

ハードウェアの違いがあるので、ファームも更新されたらしい。

そういや、こうしたコードを挙げておきながら、HP Prime実機への転送方式、気にした事がありませんでした。どなたか、ご教示ただけると幸いです。やはり、HP Conn-Kitしかないのかなァ ?

С Новым годом 2020 !

昨年もお世話になりました。今年も宜しくお願い申し上げます。

今年はロシア語で「謹賀新年」です。web時代で、簡単に見つかりますネ。

昨年末の慌ただしい中、性懲りもなくこんな計算をしておりました。

cf. ぽてん (@potenten87) 様 tweet

数学ができないとWi-Fiに繋げないホテルがロンドンにあるって。

https://twitter.com/potenten87/status/1190161683638960130

〽ロンドン、ロンドーン ! 楽しいロンドン愉快なロンドン、ロンドンロンドーン !!
(歳がバレますなァ)

この式、計算すると円周率piになるそうで、電卓でもって早速やってみた次第。

Numworks web simlatorの結果

HP Prime β版simulator

数値積分をやってみると、πに近い値になるようです。

Numworksのweb sim (https://www.numworks.com/simulator/) はJavascriptで実現されているそうで、もしかすると実機とは異なる答を出すかも知れません。

また、HP Primeのシミュレータはβ版なので、結果はこんな値でした。
XCASが使えるのですが、やり方、よく判らんですよ。申し訳ない。

量販店に行き、店頭のfx-CG50でも試してきました。
すると、CG50ではズバリπを返しますネ ! これはオモシロイ。

ただ、角度モードがRADになっている場合に限るようです。

それは、数値積分を愚直にやっていて、RADモードにした際の計算誤差がπに近いからなのだろうか、と想像します。CASIOの数値表記機能がここで威力を発揮したのでしょう。

計算の結果がπになる、ちうのは、上記のTweetをみていくと、色々出ております。奇関数、偶関数んて、すっかり忘れておりましたヨ。オツムの足りない当方なので、久しぶりに筆算をして、ようやく納得がいった次第です。今時の若い人も、こんな計算をして大学受験に備えておるのですが、大学出てしまうと、すっかり忘れてしまう人も多いかも ? (いや、ソレは当方だけなのか ?)

 

HP Primeシミュレータのβ版

昨年の2018/12/03、CASIO fx-CG50にpython機能の追加を含むファームウェアの更新情報がアップされました。
それから1年が経とうとしており、そろそろ、CASIOから新しい話題がありそうな時節となってきております。

一方、MoHPCに、HP Prime シミュレータのβ版が登場しているという話題があり、早速試してみました。

cf. HP Prime Beta Software
https://www.hpcalc.org/prime/beta/

βなので「無保証」ですが、最新のファームに対応しているらしく、3Dグラフも描けます。
ウーン、やっぱ、イイなぁ。実機が欲しくなっちゃったヨ !

先日、HP Prime G2がAmazon JPで40,000円を割る価格で出ておりましたが、既に業者も売り切ってしまったらしく、価格が元に戻ってしまいました。

最近の当方、グラフ電卓向けのコードを作ったりしておらず、もっぱら、グラフを描いて遊んだりばかりです。
こうした用途ならば、メモリが多いG2でなくてもいいのかな、などと思う所です。

HP Prime G1ですが、Amazon JPではスリーゲートがかなり安く扱っております。当方としては、スリーゲート推しでいきたい。何故か ? Amazonのページを見ると、なにやら「HP PPL超入門」というものを添付しておるのですネ。スリーゲート謹製なのかしらん、大変気になります。

https://www.amazon.co.jp/gp/offer-listing/B00OX59OYG/ref=dp_olp_new_mbc?ie=UTF8&condition=new&m=A34V9I8YOA1OXY

スリーゲートの輸入電卓には、スリーゲート謹製(?)の簡易マニュアル、読本が添付されているらしい。fx-CG50は、CASIOサイトから日本語マニュアルがダウンロードできるので、スリーゲートとしては扱う理由がなくなってしまったのか ?

スリーゲート謹製のマニュアル類、そこに、スリーゲートの電卓製品に向けられた愛を感じずにはいられないのです。

この記事を読まれた方で、スリーゲート謹製マニュアルをお持ちになられた方がいらっしゃるならば、レビューをお寄せ戴きたく !

HP Prime G2が安くなる ?

先日、HP Prime G2を4万円台で販売するショップがありましたが、今はすでに売り切れております。

HP Prime G2グラフ計算機 - Ama JP
https://www.amazon.co.jp/gp/offer-listing/B07HF6RXGG/ref=dp_olp_new?ie=UTF8&condition=new

先日、米中の貿易交渉で、一旦「和解」となったためらしく、アメリカでもHP Prime G2がおおっぴらに流通を始めた様子です。

Prime G2 available in USA - MoHPC
https://www.hpmuseum.org/forum/thread-13231.html

BestBuyでは、$150らしい。ウーム、安い !

アメリカで販売が始まったので、並行輸入業者さんも「安心価格」で流通して呉れることを強く希望します !

HP Primeでpythonが動くとな !

【重要】

この記事の内容、今日ではいささか古くなりにけり。

以下のblog様の記事が参考になります。

 

ref. HP PrimeのPython覚え書きなのだ！
http://hangyodon.cocolog-nifty.com/hangyodon/2022/01/post-374e72.html
 

  ---

 

先日、Sentaro 様からお知らせ戴きまして、HP PrimeのXCASがpython文法に対応した事で、Primeでもpythonコードが動く様になったとの事。手元で動いているHP Primeシミュレータでpythonコードを動かしてみましたヨ。

以下、少々長いのですが。programエディタで「CAS」のチェックボックスをONにしてから入力して下さい。

#cas

## title : py_shida.py
## begin : 2019-02-23 18:44:05 

def py_shida() :
  ## window size
  WIDTH  = 320
  HEIGHT = 240

  ##  field 
  XL  = -4.0
  XH  =  4.0
  YL  =  0.0
  YH  = 11.0

  ##  反復計算の回数
  count  = 100000

  ##  シダの葉を描くよん
  cnt  = 0
  r    = 0
  px   = 0
  py   = 0
  
  while cnt < count :
    cnt = cnt + 1
    r = RANDOM()
    if   (r < 0.01) :
      x = 0
      y = 0.16*py
    elif (r >= 0.01 and r < 0.07) :
      x = 0.2*px  - 0.26*py
      y = 0.23*px + 0.22*py + 1.6
    elif (r >= 0.07 and r < 0.15) :
      x = -0.15*px + 0.28*py
      y =  0.26*px + 0.24*py + 0.44
    else :
      x =   0.85*px + 0.04*py
      y =  -0.04*px + 0.85*py + 1.6
  
    px = x
    py = y
    i = (py-YL)/(YH-YL)*WIDTH
    j = HEIGHT-(XH-px)/(XH-XL)*HEIGHT
    PIXON_P(i, j, RGB(0, 96, 0))
 

  ##  key wait
  FREEZE

#end

ンマァ、なんて長いのでしょう ! まるで某代議士のチ○ポの様な ... オゲレツ !!

それはさておきまして、OSを最新にすると、XCASがpython文法に対応したとの事で、HP Primeでpythonコードが動く様になります。HP Primeをお持ちの方は、是非ともお試し頂きたく !

なお、βではありますが、新しいシミュレータなどが、つぎの所から入手できます。


HP Prime Beta Software - MoHPC
https://www.hpcalc.org/prime/beta/


一方、海の向こうではGeneration 2nd という新型ハードウェアのHP Primeが販売されているとか。性能が向上していると ! 欲しくなっちゃったヨウ。

で、並行輸入品がないのか、調べてみましたら、一応はあるらしい。

HP Prime G2グラフ計算機 - Amazon JP
https://www.amazon.co.jp/%E3%83%92%E3%83%A5%E3%83%BC%E3%83%AC%E3%83%83%E3%83%88%E3%83%BB%E3%83%91%E3%83%83%E3%82%AB%E3%83%BC%E3%83%89-2AP18AA-B1S-Prime-G2%E3%82%B0%E3%83%A9%E3%83%95%E8%A8%88%E7%AE%97%E6%A9%9F/dp/B07HF6RXGG/ref=sr_1_6?ie=UTF8&qid=1552733648&sr=8-6&keywords=hp+prime+graphing+calculator

価格: ￥ 64,920 ... 高いなぁ。
しかも、アヤシイのが「256 MBフラッシュメモリ」という所。確か、G2はFlash 512 MBだったような ?

HP Prime - en Wikipedia
https://en.wikipedia.org/wiki/HP_Prime

それでも「2AP18AA#B1S」という型番が出ているので、つられてしまうのです。さて、その真相は !?

2019 Prosit Neujahr !

新年・2019年を迎えました。遅ればせながら、謹賀新年。

2018年は、高機能電卓界隈でも動きがありました。

1. CASIO fx-CG50で、MicroPythonが動作する様になりました
   
   当初、フランスの教育関係者向けに限定供給されていた、GRAPH 90+E 向けのMicroPythonがfx-CG50向けにも提供、OSのアップグレードで利用できる様になりました。
   
   
2. HP Prime のハードウェアが更新される
   
   Sentaro様にお知らせ戴いたのですが、HP Primeの第二世代ハードウェアが発売されている様です。
   
   cf. New G2 HP Prime - MoHPC
   http://www.hpmuseum.org/forum/thread-11202.html
   
   どうやら、MPUのクロックが高くなり、処理速度が向上している模様。
   しかしながら、流通がこれかららしく、米国でも入手が容易ではない状況です。
   
   また、最新のファームで、CAS環境でMicroPythonが使えるとの事。これは新旧、両方のハードウェアで利用可能の模様。
   
3. フランスNumworks社のMicroPython電卓
   
   既に、2016年の時点で、海外では仏Numworks社のグラフ電卓が発売されているとの事。
   
   cf. NumWorks - The graphing calculator that makes learning math easier
   https://www.numworks.com/
   
   Numworksの特徴は、MicroPythonが動作する所でしょうか。計算機能などもかなりユニークです。グラフ電卓なので、関数の2D plotや、統計グラフ機能などがあります。複素数や行列の計算も用意されております。
   電卓の機能としてはかなりのものですが、各アプリケーション間でデータの共有ができないらしく、惜しい印象です。
   一方、ハードウェアの設計、OSを含めたソフトウェアも、全て公開されているので、ユーザーが独自の拡張を行う事が、原理的には可能らしい。とても興味深い製品ではありますが、果たして(高機能)「電卓」なのか ? という疑問もあります。
   オープンソースで、サイトにはブラウザで動作するシミュレータも用意されてます。興味のある方は、お試しいただきたく。
   

2019年も始まったばかりで、これから何があるのか。不安でもあり、楽しみでもあります。

今年も宜しくお願い申し上げます。

「シダの葉グラフィクス」

 2017-10-20 には、CASIO fx-CG50が国内発売される、との事で、カラーグラフィクスの表示できる高機能電卓向けに、「シダの葉グラフィクス」のプログラムを実行してみました。

シダの葉グラフィクスについては、つぎを御参照。

反復関数系 - wikipedia
https://ja.wikipedia.org/wiki/%E5%8F%8D%E5%BE%A9%E9%96%A2%E6%95%B0%E7%B3%BB


当方、未だカラーグラフィクスの高機能電卓を持っていないので、HP Primeのシミュレータで動作するコードを作成してみました。
以下の版は、モノクログラフィクスでの表示ですが、RGB()のパラメタを少し変えるだけで、カラーになります。

プログラム「draw_shida.ppl」

//  title : draw_shida.ppl
//  begin : 2017-10-03 15:18:44 
//  

//  
width  := 320;
height := 240;

//  initial ooord
x := 0;
y := 0;
//  loop counte
cnt :=  30000;
//  working
tmp := 0;

xmin := -3;
xmax :=  5;
ymin :=  0;
ymax :=  8;

//  functions

//  f1()
f1()
BEGIN
  tmp := 0;
  y := 0.16*y;
  x := tmp;
END;

//  f1()
f2()
BEGIN
  tmp := 0.2*x  - 0.26*y;
  y := 0.23*x + 0.22*y + 1.6;
  x := tmp;
END;

//  f1()
f3()
BEGIN
  tmp := -0.15*x + 0.28*y;
  y :=  0.26*x + 0.24*y + 0.44;
  x := tmp;
END;

//  f1()
f4()
BEGIN
  tmp :=  0.85*x + 0.04*y;
  y :=  -0.04*x + 0.85*y + 1.6;
  x := tmp;
END;


//  MAIN 
EXPORT DRAW_SHIDA()
BEGIN
// Clean the screen (G0):
RECT();

  LOCAL  i, r, col;
  LOCAL  xp, yp;

  FOR i FROM 1 TO cnt DO
    r := random();
    IF (r < 0.01)               THEN f1(); col := RGB(0, 0, 0); END;
    IF (r >= 0.01 AND r < 0.07) THEN f2(); col := RGB(0, 0, 0); END;
    IF (r >= 0.07 AND r < 0.15) THEN f3(); col := RGB(0, 0, 0); END;
    IF (r >= 0.15)              THEN f4(); col := RGB(0, 0, 0); END;
    xp := (x-xmin)/(xmax-xmin)*width;
    yp := (y-ymin)/(ymax-ymin)*height;
    PIXON_P(yp, xp, col);
  END;


//  key wait loop
REPEAT UNTIL GETKEY() == -1;
FREEZE;

END; 

小汚いコードですが、 移植は簡単だと思います。

懸垂線、再び

以前、当blogで、関数電卓に付き物の双曲線関数の使用例で、懸垂線の計算について調べたネタをアップしました。

懸垂線は、紐や鎖の両端を持ってタラリと垂らすと、その紐なり鎖なりが描く曲線なんだそうです。数式では、つぎの式で示されます。

y = a * cosh(x/a)

HP35SのTraining modules には、懸垂線の計算例が紹介されておりましたが、紐の両端を持つ、その高さが等しい例でした。
しかし、紐をタラリと垂らすとなれば、その両端の高度が違う場合もありましょうや。
そこで両端の高度が異なる場合、懸垂線のパラメタ (カテナリ数と極点x座標)の数値を求める方法を調べてみようと思います。

懸垂線は、紐の両端の高度が異なっても、同じ式で表現されると言います。
両端の高度が異なる懸垂線は、煎じ詰めるとつぎの数値で決定される、として良いでしょう。

• 紐の片方の端を原点とし、もう片方の端を「吊り下げ点」として、吊り下げ点の座標(D, H)を与える
• また、紐の長さ L を指定する

こうすれば、懸垂線のパラメタは一意に決定されそうです。これら3つのパラメタ、H, D, L を指定すれば、懸垂線のパラメタである a (カテナリ数), x (極点x座標)が決定され、懸垂線自体が決定できるだろうというハラです。
コチャコチャと計算すると、D, H, Lからa, xを決定する方程式は、つぎの2本である事が判ります。

a*(sinh(x/a)+sinh((D-x)/a)-L = 0
a*(cosh((D-x)/a)-cosh(x/a))=H = 0

非線形の連立方程式です。解析的に解けるんかいや ?
大変難しそうで、当方の知見は及びません。しかし、高機能電卓の一部には、数値解を計算する機能を持ち合わせているものがあります。そう、HP PrimeやHP50Gが、まさにソレなのです ! (もちろん他にもありますが、当blogはHP電卓の情報発信を任じておりますので)。

HP50Gでは、つぎの様にします。

1. D, H, L に数値をストアします
2. つぎの通り、スタックに配します



[ 'A*(SINH(X/A)+SINH((D-X)/A)-L=0' 'A*(COSH((D-X)/A)-COSH(X/A))=H=0' ]
[ 'A' 'X' ]
[ 10 10 ]

3. MSLVを実行

最後にスタックに積んだ [ 10 10 ] は、数値解を求めるための「初期値」で、未知数A, X に対応します。この例ではA=!0, X=10 としましたが、D, H, Lとの兼ね合いから、適切な値を指定しないと数値解が得られない事もあります。

例として、つぎの場合を計算させてみましょう。

入力  D=35, H=18, L=44

結果  A=19.0231959735, X=9.23406136377

HP Primeの場合は、もっと楽 ?

「解く」アプリを呼び出し、つぎの様に数式を入れます。

そして、[Num」キーを押して、D, H, Lに数値を設定、チェックボックスのチェックを外し、A, Xを求めるので、A, Xのチェックボックスにはチェックを入れたままにして、スクリーンキー右端の「解く」を押します。

答え。

エミュレータでは、瞬時に値を得られましたが、実機では若干、時間がかかるでしょうか。それでも、50Gよりは断然早いものと思われます。

電卓による母比率計算の例

先般、参議院議員選挙が行われ、改憲勢力が2/3を獲得したとの事で「平和憲法サヨウナラ」となる事が決まりました。

マスコミはこぞって「出口調査の結果に拠りますと」と述べておりますが、統計学の教える所に拠りますと、偏りのない標本集団を得られれば、母集団の比率が推測できるのだと言います。

HP50GやHP Prime、TI-83+には、母比率の推定機能が盛り込まれているので、手のひらで計算する事が出来ます。

例えば、標本数n=320の内、あるカテゴリの数がx=170だった、として、母集団における、あるカテゴリの割合=母比率がどのくらいになるのか、計算できるのです。

標本集団における、あるカテゴリの比率はx/n = 0.53125 ですが、これは標本集団の比率であります。

標本集団の標本数を十分大きく取る事で、標本比率は正規分布に近付いていく、という難しい理屈(「中心極限定理」)があるのだそうですが、電卓の機能は、そうした背景を以って構成されているので、単純に数値を入力してやれば、立ちどころに答が得られます。

しかし、背景となる理屈の詳細までは知らずとも、ある程度の知識を持っていないと、その答えを吟味する事は難しいものです。

母集団から一部の標本を抽出し標本比率を計算すると、その標本比率は、やはり正規分布をする、のだそうです。

標本比率の分布する範囲は、正規分布となるため、例えば95%程度の確度で、この範囲に存在する、という結果が得られる事になります。

n=320, x=170の例で、確度=0.95 (95%)の分布範囲は、0.4765744 〜 0.5859256 の範囲、という結果が返ってきます。

実際に、HP Primeでやってみましょう。

[Apps] -> [推論] とキーを押して、つぎの情報を選択します。

• 方式 : 信頼区間
• タイプ: Z-Int: 1 π

[Num]キー押下、次の数値を入力します。

• x:170
• n:320
• C:0.95

これで計算の準備が出来ました。つぎの選択肢があります。

• [計算]メニューキーを押下すると、結果ウィンドウで表計算っぽい結果が得られます
• [Plot]キーを押下すると、正規分布グラフと共に結果が表示されます

Plotキーを押してグラフを表示した例

HP Prime Virtual calcの新版があるらしいが

ftp://ftp.hp.com/pub/calculators/Prime/

ここを見ると、新しいタイムスタンプのファイルがありますが、ファイル名が「去年」です。
しかし、念の為、アップデートを行ってみました。先ずは、仮想電卓を動かして更新を掛けますが ... ダメ。
そこで、ファイルをダウンロードして実行してみますが ... 「新しいバージョンがインストール済み」と出て、失敗します。
ウーン、どうなっておるのか !?

どうやら、日付を間違えて設定されているのではなかろうか ? 近々、修正されてアップされるのではないか、と密かに期待します。

速報・HP Prime仮想電卓の新バージョン

New Version: 2016.04.14 r10077 - MoHPC
http://www.hpmuseum.org/forum/thread-6085.html

という訳で、仮想電卓の新版も出ている様です。
しかし、仮想電卓を動かしても、アップデートされる気配がありません。ftpのファイル置き場を覗いたら、どうも新旧2つのファイルがあり、アップデータが混乱しているのか ?

ftp://ftp.hp.com/pub/calculators/Prime/

HP Primeによる複素数のグラフの例

最近はJulyさんのAmazon店で、HP50Gの在庫が5台と少なくなりつつありますが、HP Primeの方は在庫も潤沢な様です。
そこで、Primeの売上に貢献すべく(?)、こんな話題をお届け。
 
HP50G, HP Primeでも複素数の計算が出来ます。

ここでは、HP Primeにて等角写像の例として、Zhukovski 変換を使ったグラフィクスを描いてみましょう。
詳しい内容は、Zhukovski変換を調べてください。

先ずは、変換する前の円の中心座標を定めます。複素平面に描くため、円の中心座標も複素数として設定しますから、複素数を保持できる変数が必要です。デフォルトで利用できる複素数変数にはZ0～Z9までがあり、ここではZ0に(-0.08, 0.15)を設定してみました。
つぎに、変換する円の半径ですが、これはsqrt((1-Re(Z0))^2+(Im(Z0))^2) とします。

変換する円の式をdefineによって定義しておきます。

更に、Zhukovski変換の式もdefineで定義すると、楽に進みます。

最後は、「パラメトリック・プロット」アプリケーションを呼び出し、つぎの様にグラフの式を設定しておきます。アプリケーションの複製を作り、Zhukovskiとでも名付けてやってもいいかも知れません。

後は、Graphキーを押してやれば、変換前の円と、Zhukovski変換にて得られる「翼断面」の図が得られるはずです。

同様の事をHP50Gでも作業できますが、実行速度はHP Primeの方が断然早そうです。