Scilab/MATLAB用xorshift関数（擬似乱数生成関数）

　Scilab/MATLAB用xorshift関数（擬似乱数生成関数）

| ← back |

| home |

前置き

Scilab/MATLABは大変有用なシミュレーション・ツールで、Scilab/MATLABでアルゴリズムを開発してからC言語に移植することもあれば、C言語のプログラムの動作検証のためにScilab/MATLABのプログラムを作成して処理結果を比較検討することもあります．
このような用途のために、Scilab/MATLABとC言語で同一の系列を出力可能、かつ任意の初期値（seed）を設定可能なxorshift擬似乱数生成関数を作成しました．
C言語で作成したxorshift関数のDLLをScilab/MATLABにリンクして使えれば良いのですが、Windows版のScilabではうまく出来なかったので少し工夫してオリジナルのScilab用xorshiftを作った次第です．
Scilabで動作確認をおこなっていますが、MATLABでも基本的には問題なく動くはずです．（ただし一部修正は必要です）　Octaveへも移植も出来ると思います．
Scilab/MATLABはインタープリタなので、動作速度は非常に遅いです．

Scilab用xorshift

Scilab等のシミュレーション・ソフトは倍精度演算を基本としていて、整数型の変数やビットごとの排他的論理和演算子（xor/exor）は無いようです．　そこで、苦肉の策として32bit整数を値-1/1を持つデータ長32の配列で表現することにしました．　例えば16進表記の 0xA5A50F0F は

　[1 -1 1 -1 -1 1 -1 1 1 -1 1 -1 -1 1 -1 1 -1 -1 -1 -1 1 1 1 1 -1 -1 -1 -1 1 1 1 1 ]

として扱います．　同様の方法で64bit/128bit整数も表現可能です．（64bit/128bitのxorshiftも作れます）
0/1の代わりに-1/1の2値表現を用いているため、排他的論理和の演算は乗算と符号反転（減算）であらわされます．（下図参照）
ソフトウェアの世界では2値信号は論理値0/1（True/False）ですが、工学分野ではこのように-1/1の2値（整数）として扱うことも多いです．
Scilab用xorshift　sci_xorshift.sce
内部状態（レジスタ）を保持するために　[y, seed]=xorshift32bit(seed) の形で関数呼び出しをおこなってください．（y：生成された乱数、seed：レジスタ／初期値）
詳しいコメントはついていませんが、コードを読めば簡単に処理内容は理解できるはずです．
関数xorshift()はMATLAB, Octaveでもそのまま動作すると思いますが、初期値（seed）計算に使う関数hex2seed()は修正が必要かもしれません．　Scilabの関数ascii()に相当する処理は、MATLABでは関数double()で出来るようです．　他にmprintf()の関数名変更も必要です．

/* Scilabe 32bit xorshift function */
/* (c) 2025 cepstrum.co.jp         */

clear;

DATLEN =80000;
HEXSEED='A5A50F0F';

//---------------------------------
// 32bit xorshift
//   p^=(p<<5); p^=(p>>13); p^=(p<<6);
//   useage : [y, seed]=xorshift32bit(seed)
//   seed is (1:32) array of -1/1 instaed of 0/1 value

function [y, seed]=xorshift32bit(x)
  x=-x.*[x(6:32), -ones(1:5)]
  x=-x.*[-ones(1:13), x(1:19)]
  x=-x.*[x(7:32), -ones(1:6)]
  seed=x
  y=sum(2^(31:-1:0).*((x+1)/2));   // convert (1:32) array of -1/1 value to unsigned integer
endfunction

//---------------------------------
// covert hex string to (1:32) array of -1/1 value

function y=hex2seed(x)
  x=hex2dec(x);
  x=dec2bin(x);
  x=ascii(x);
  x=2*(x-48)-1;      // 48 (0x30) is ASCII code of '0'
  y=x
endfunction

//---------------------------------

seed=hex2seed(HEXSEED);    // generate xorshift seed

mprintf("seed(uint) %u\n", hex2dec(HEXSEED));
mprintf("seed(-1/1) ");
for i=1:32
  mprintf("%2i ", seed(i));
end
mprintf("\n");

y=zeros(1:DATLEN);
for i=1:DATLEN
  [y(i), seed]=xorshift32bit(seed);    // seed holds internal state
end

for i=0:7
  mprintf('%12u %12u %12u %12u\n', y(i*4+1:i*4+4));
end

savewave('out.wav', y/max(abs(y))-0.5, 8000);

実行結果

seed(uint) 2779057935
seed(-1/1)  1 -1  1 -1 -1  1 -1  1  1 -1  1 -1 -1  1 -1  1 -1 -1 -1 -1  1  1  1  1 -1 -1 -1 -1  1  1  1  1 
  1344132296   3134319730   3660722449   2152518865
  4250013827   2478612643   4023055366   3013925157
   155769594    888360614    225206834   2578317235
  2848799107   2345412211   1636953940   3087578548
  4293475012   2975037943   1428104117   1237602575
   573791644   1923676699   3390551366   2981361934
  1978748072   1014979822   1472610184   1848186630
  1755662809   2146120730    551644919    261671666

C用xorshift

C用xorshift　c_xorshift.c
Scilab用のプログラム sci_xorshift.sce と同一の擬似乱数が生成されます．

/* 32bit xorshift function */
/* (c) 2025 cepstrum.co.jp */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdint.h>

#define SEED      0xA5A50F0F
#define DATLEN    80000
#define OUT_FNAME "out.txt"  // output data file name

//================================================================
// xorshift random number generator (32bit)

uint32_t xorshift32bit(uint32_t seed) {
  static uint32_t y=0x0f0f0f0f;    // default seed value

  if (seed!=0) {
    y=seed;                        // set new seed
    return seed;
  }

  y=y^(y<<5);
  y=y^(y>>13);
  y=y^(y<<6);
  return y;
}

//================================================================

unsigned rnd[DATLEN];

int main() {
  FILE     *fp;
  unsigned i, k;

  xorshift32bit(SEED);     // set seed

  printf("seed %u\n", SEED);

  for (i=0; i<DATLEN; i=i+1) rnd[i]=xorshift32bit(0);

  for (i=0; i<8; i=i+1) {
    for (k=0; k<4; k=k+1) {
      printf("%12u ", rnd[i*4+k]);
    }
    printf("\n");
  }

  // output xorshift random noise
  fp=fopen(OUT_FNAME, "w");
  for (i=0; i<DATLEN; i=i+1) fprintf(fp, "%u\n", rnd[i]);
  fclose(fp);

}

実行結果

seed 2779057935
  1344132296   3134319730   3660722449   2152518865 
  4250013827   2478612643   4023055366   3013925157 
   155769594    888360614    225206834   2578317235 
  2848799107   2345412211   1636953940   3087578548 
  4293475012   2975037943   1428104117   1237602575 
   573791644   1923676699   3390551366   2981361934 
  1978748072   1014979822   1472610184   1848186630 
  1755662809   2146120730    551644919    261671666

| ← back | ↑top |

| home |