Sinテーブルを極限まで（？）小さくする。

なんとなくtiny13でどーにかしようとしているドM野郎の冨塚です。

要は「マルツtiny85だと秋月tiny13が４つ買える。」ってことです。
（メモリは８倍なので得なのだけど）

全然電子工作エリアではないのですが
　如何にゼニをかけないかという点で
　　セコセコと最小構成で遊びたいわけです。

○三角関数系でいかにメモリを小さくするかで色々考えてたら思い出しました。

以下のネタは、1980年代の8bitパソコンで３次元ワイヤーフレームアニメーションをやるため、いかに小さなメモリで三角関数を計算するか考えてマシン語で組んでたのがルーツです。
（なんせFM-7のサブシステム共有メモリは128バイトしかなかったのです、RGBの一枚を捨てて、VRAM上でプログラムを走らせてました。）

☆テーブルは60バイトでなんとかなっちまう。
　60バイトってのが重要で、tiny13のEEPROMエリアが64バイト。
　つまりフラッシュメモリの1kbyteをプログラム領域としてフルに使用できるのね。

☆0度から360度の１度毎のdegで精度はほぼ16bitです。

○ 下準備
１、0度から90度までのSin値を32768倍したものを10度毎に
　　20バイトのメモリ上に置きます。
　　（なぜこいつだけ32768にするのかは後述）

２、０度から9度までのSin値を65536倍したものを1度毎に
　　20バイトのメモリ上に置きます

３、０度から9度までのCos値を65536倍したものを1度毎に
　　20バイトのメモリ上に置きます

C言語的に書くとこんなね。数値はホンマにつかえます。

uint16_t sin90[10] = { 0,5690,11207,16383,21062,25101,28377,30791,32270,32768};
uint16_t sin10[10] = { 0,1143,2287,3429,4571,5711,6850,7986,9120,10252};
uint16_t cos10[10] = { 0,65526,65496,65446,65376,65286,65176,65047,64898,64729};

 このunsigned16bitの30配列だけで0-360度の三角関数をなんとかします。
まずSin値を求めたい角度の１の位を捨ててφとし、1の位をθとします。

４、先ずSinφとCosφを求めます。
○10度毎のSinφテーブルは90度までしかないのですが
100-180度はSin(180-φ)　190-270度は-Sin(φ-180)　280-350度は-Sin(360-φ)
で求めます。

○10度毎のCosφは上記に90度足してあげれば求まります。

５、次に１の位のSinθとCosθを求めます。
これはそれぞれ０度から９度までしかないので、テーブルを参照すればOk

６、この公式を使います。
     Sin(φ+θ) = Sinφ×Cosθ+Cosφ×Sinθ
　これで０度から360度のSin値が求まるわけです。
　
☆なぜ 10度毎のSinφの値が32768倍なのか。
　つまるところ整数値なので
　(Sinφ×65536)×(Cosθ×65536)+(Cosφ×65536)×(Sinθ×65536) / 65536;
   でいいのですが、これやっちゃうと計算中に32ビットで桁溢れします。
　まあlong64でやればいいんですけど、8bitCPUで64ビット演算なんかあんまりしたくないしバイナリコードがでかくなったら本末転倒。
　精度はオチるがSinφのテーブルは32768倍にして
　(Sinφ×32768)×(Cosθ×65536)+(Cosφ×32768)×(Sinθ×65536) / 32768;
　とやることでint32で計算できるというわけです。

2π=360でなく2π=256でやるなら
上位5bit（16バイト）下位3bit（32バイト）計48バイトでいけるか?

サンプルはgccでコンパイルできる普通のcのソースです。
※tiny13用ではありません。

#include <stdio.h>
#include <math.h>
#define PI 3.141592

unsigned short  sin90[10] = {  0, 5690, 11207, 16383, 21062, 25101, 28377, 30791, 32270, 32768};
unsigned short  sin10[10] = { 0,1143,2287,3429,4571,5711,6850,7986,9120,10252};
unsigned short  cos10[10] = { 0,65526,65496,65446,65376,65286,65176,65047,64898,64729};

int sin36(int d ) {
    d = d % 36;
    if ( d <= 9 ) return sin90[d];
    if ( d <= 18 ) return sin90[18-d];
    if ( d <= 27 ) return -sin90[d-18];
    if ( d <= 36 ) return -sin90[36-d];
}

int cos36(int d) {
    return sin36(d+9);
}

int sin16(int d) {
    int d36,d10,value ;
    int s36,s010,c36,c010;
    d = d % 360;
    d10 = d % 10;
    d36 = ( d - d10 ) / 10 ;
    s36 = sin36(d36) ;
    if ( d10 == 0 ) return s36*2;
    c36 = cos36(d36) ;
    s010 = (int)sin10[d10] ;
    c010 = (int)cos10[d10] ;
    value = (s36 * c010 + c36 * s010) /32768 ;
    return (int) value ;
}

int cos16(int d) {
    return sin16(d+90);
}

int main(void)
{
    int d , s16 , c16;
    double deg , rad , rsin , tsin, rcos ,tcos;
    double dsin ,dcos,dmaxs,dmaxc ;
   
    /*どのくらいの精度が出るのか、math.hと比較*/
    dmaxs = dmaxc = 0.0;
    for ( d = 0 ; d < 360 ; d+=1 ) {
        deg = (double) d;
          rad = deg * PI / 180.0;
          rsin = sin(rad);
        s16 = sin16(d) ;
        tsin = (double)s16 / 65536.0 ;
          rcos = cos(rad);
        c16 = cos16(d) ;
        tcos = (double)c16 / 65536.0 ;
        dsin = rsin - tsin ;
        if ( dsin < 0 ) dsin = -dsin ;
        if ( dmaxs < dsin ) dmaxs = dsin ;
        dcos = rcos - tcos ;
        if ( dcos < 0 ) dcos = -dcos ;
        if ( dmaxc < dcos ) dmaxc = dcos ;
        printf("sin(%.3d)  math.h[%1.6f]table[%1.6f]Gosa[%1.6f] cos(%.3d) math.h[%1.6f]table[%1.6f]Gosa[%1.6f]\n"
        ,d,rsin,tsin,dsin,d,rcos,tcos,dcos);
    }
    printf("max Gosa sin:%1.6f  cos:%1.6f\n",dmaxs,dmaxc);
      return 0;
}