sin of calc

Q: /usr/include/mth.h に パイの定数が M_PI として定義されてます。 3.14159265358979323846 これ小数点を除くと21桁です。 一方、ieee754の定義では、仮数部22Bit で、10進数換算16桁となります。 ソースの過剰定義と思えてなりません。何か理由が有るのですか？

A: `<mth.h>` における `M_PI` の定義が小数点以下の桁数で IEEE 754 に準拠していないのは、確かに興味深い点ですね。一般的には、数学定数としての円周率 π は 3.14159265358979323846 と表されますが、IEEE 754 の倍精度浮動小数点数では約 15 〜 16 桁の精度しか保証されません。

一部のライブラリやヘッダーファイルでは、これを補完的な定数として提供していることがありますが、これは厳密な浮動小数点数の精度とは異なります。

このような過剰な桁数を定義する理由は、主に数学的な精度や扱いやすさを考慮した結果として考えられます。C 言語の `<mth.h>` など、特定のヘッダーファイルやライブラリでは、数学的な定数を便宜上拡張して提供することがあります。これは、プログラマがその精度を要求する場合には適切に丸めることが期待されるものです。

また、特定のプラットフォームやコンパイラに依存する場合もあります。そのプラットフォームやコンパイラが、自身のニーズに合わせて定義を行っている可能性も考えられます。

要するに、C 言語のヘッダーファイル内の定数は、数学的な定義とプログラミング上の実装のバランスを取ったものであり、必ずしも厳密な浮動小数点数の規格に合致する必要はない場合があります。

doubleの有効桁数が16桁と知ったけど、それってドンダケー凄いの？ 身の周 りで実感してみる。そう、世界が100人の村だったらって本が昔流行したみた いに、身辺感覚で捉えてみるのさ。

今の時代、人生100年と良く言われる。生まれた時に、蝋燭に点灯。それが燃 えつきたら、天寿を全うしたって奴ね。その蝋燭には、秒の刻みがついている。 さて、その刻みは、いか程に。また、何人分の刻みが、16桁で実現できる？

schemeでササッと計算。

gosh$ (/. (expt 2 52) (* 3600 24 365 100))
1428082.0736207813

goshを入れていない人向けに、bc で、検算する例をば。ドットは直前の演算 結果を保持しています。

3600 * 24 * 365 * 100
3153600000
2 ^ 52 / .
1428082
. / 10000
142

142万人分の人生を刻めるとな。とんでもない桁数だ。2021年の100才以上の日 本の人 口は、86510人らしい。また、都市人口だと、ピッタリは流石にない。京都市 が146万人で一番近いかな。次点はさいたま市の132万人あたりか。

今迄、余り気にしなかったけど、パッケージの検索順ってどうなってるのだろ う？

 *** Welcome to IELM ***  Type (describe-mode) or press C-h m for help.
ELISP> load-path
("~/.emacs.d/lisp" "/home/sakae/.emacs.d/elpa/cargo-20230125.1253" ...

ちょっと醜いので、ielmな結果をファイルに書き出して整理してみた。

[sakae@fb /tmp]$ cat path | tr ' ' '\n'
("~/.emacs.d/lisp"
"/home/sakae/.emacs.d/elpa/cargo-20230125.1253"
"/home/sakae/.emacs.d/elpa/cargo-mode-20230213.401"
  :
"/usr/local/share/emacs/29.1/lisp/calendar"
"/usr/local/share/emacs/29.1/lisp/calc"
"/usr/local/share/emacs/29.1/lisp/obsolete")

ほー、calcは最後に検索されるっていう不遇な位置になってるのね。

こんな風にして、配布された物のコピーを作る。

vbox$ mkdir calc
vbox$ cd calc/
vbox$ cp /usr/local/share/emacs/29.1/lisp/calc/*.el.gz .
vbox$ for f in *
> do
> gunzip $f
> done

忘れないように下記を、.emacs/init.el に追加。

(setq load-path (append '("~/calc" ) load-path))

これで検索順が一番になるはず。

(defun calcFunc-sin (x)   ; [N N] [Public]
=>(cond ((and (integerp x)
              (if (eq calc-angle-mode 'deg)
                  (= (% x 90) 0)
                (= x 0)))

おお、来た来た。丁度 x を評価し終ったんだけど、mini-bufferには、 Result: 45 (#o55, #x2d, ?-) と、入力した値が表示されてる。

歩を進めると、一瞬で、

        (t (calc-record-why 'scalarp x)
=>         (list 'calcFunc-sin x))))

に到達して、mini-bufferには、こんな結果が表示された。

Result: (float 707106781187 -12)

もういちど、のろまに進んでみるか。 それよりもキビキビとソースの海を泳いでみるか。

TAGSを作成すれば、どんどん関数間を移動できる。使わない手はないな。

vbox$ etags --help
  :
  c          .c .h
  c++        .C .c++ .cc .cpp .cxx .H .h++ .hh .hpp .hxx .M .pdb
  c*         .cs .hs
  cobol      .COB .cob
  lisp       .cl .clisp .el .l .lisp .LSP .lsp .ml
   :
  python     .py
  ruby       Rakefile Thorfile .rb .ru .rbw
  rust       .rs
  scheme     .oak .rkt .sch .scheme .SCM .scm .SM .sm .ss .t
   :
vbox$ etags -V
etags (GNU Emacs 29.1)
Copyright (C) 2023 Free Software Foundation, Inc.
This program is distributed under the terms in ETAGS.README

昔と違って多様な言語をサポートしているな。ctagsグループの支援を得てい るのでしょうか。

etags *.el して、すいすいとソースの海を徘徊できる様になった。あちこち をウロウロしたけど泳ぎ 疲れて、結局は核心に迫れず。

ここで一つ注意がある。load-pathに登録されてるソースは、所定の場所にあ ると想定されるTAGSを参照するって事。だから、規定の場所からソースを取り 出してきてTAGSを作成しても、それは参照されない。結果として、そんな関数 は登録されていないとなってしまう。

もう一つ注意。補完すると、こんなのが出現。これしか、飛んでいけない。

In this buffer, type RET to select the completion near point.

13 possible completions:
calc                      calc-dispatch             calc-embedded
calc-embedded-activate    calc-eval                 calc-grab-rectangle
calc-grab-region          calc-grab-sum-across      calc-grab-sum-down
calc-keypad               calc-undo                 calculate-lisp-indent
calculator

自由にあちこち飛ぶには、calcを一度実行する必要があるみたい。そうすると 多数の候補が出現して、よりどりみどりになる。

これって、パッケージがロードされると、それも対象にしてくれるからかな？ math- で始まる語句は、どうやらcalcを下支えするプリミティブみたいだけど、 これも堂々とノミネートされるよ。

折角etagsが出てきたので、原点を探ってみる。ctagsとペアね。

[sakae@deb lib-src]$ wc ctags.c etags.c
     2      5     35 ctags.c
  8071  30615 211848 etags.c
  8073  30620 211883 total
[sakae@deb lib-src]$ cat ctags.c
#define CTAGS 1
#include "etags.c"

うり二つのソーセージもとえ双生児だわい。で、中身のrubyの項でも散策。

static void Prolog_functions (FILE *);
static void Python_functions (FILE *);
static void Ruby_functions (FILE *);
static void Rust_entries (FILE *);
static void Scheme_functions (FILE *);

static const char *Ruby_suffixes [] =
  { "rb", "ru", "rbw", NULL };
static const char *Ruby_filenames [] =
  { "Rakefile", "Thorfile", NULL };
static const char Ruby_help [] =
  "In Ruby code, 'def' or 'class' or 'module' at the beginning of\n\
a line generate a tag.  Constants also generate a tag.";
static const char *Ruby_interpreters [] =
  { "ruby", NULL };

static language lang_names [] =
{
    :
  { "ruby",      Ruby_help,      Ruby_functions,    Ruby_suffixes,
                 Ruby_filenames, Ruby_interpreters },

/*
 * Ruby support
 * Original code by Xi Lu <lx@shellcodes.org> (2015)
 */
static void
Ruby_functions (FILE *inf)
{

これ、rubyで検索して出てきたもの。肝心の本体部分は、165行であった。 pythonは誰の寄贈かと調べたら、Idea by Eric S. Raymond だってさ。見てい て、あきないな。

17.5.5.2 The Sine Function に ピッタリの例が出てた。

(defmath mysin ((float (anglep x)))
       (interactive 1 "mysn")
       (setq x (to-radians x))    ; Convert from current angular mode.
       (let ((sum x)              ; Initial term of Taylor expansion of sin.
             newsum
             (nfact 1)            ; "nfact" equals "n" factorial at all times.
             (xnegsqr :"-(x^2)")) ; "xnegsqr" equals -x^2.
         (for ((n 3 100 2))       ; Upper limit of 100 is a good precaution.
           (working "mysin" sum)  ; Display "Working" message, if enabled.
           (setq nfact (* nfact (1- n) n)
                 x (* x xnegsqr)
                 newsum (+ sum (/ x nfact)))
           (if (~= newsum sum)    ; If newsum is "nearly equal to" sum,
               (break))           ;  then we are done.
           (setq sum newsum))
         sum))

ここでも、テーラー展開である。defmathの中では、for構文が使える。これで、 展開の次数をどんどんと増加しながら計算。前の結果とほぼ等しかったら(浮 動小数点なんで、等しいかの比較は不可)終了という戦略だ。

これに続く節では、色々と注意点が説明されてた。

このマクロがどう展開されるかの補助関数は、残念ながら提供されていない。

うえのコードを実際に確認できる。.emacs/calc.el ってファイルに記述す るだけ。

'mysin(45) と入力。代数式って扱かいだ。

alg' mysin(45)
     0.707106781179

calcの右窓に、こんな結果が出てくる。

別の方法として、この関数をキー・バインドしちゃって、計算したい値をスタッ クに置き、それからキー操作するって事もできる。

まずは、キーの登録。ZDと入力すると、下記の質問がなされる。

Define user key: z-q
Set key z q to command: calc-mysin

使ってみる。もう普通の気分で、30 zq するだけだ。

     30
mysn 0.499999999998