STM32VLDiscoveryで100の階乗計算をさせてみた

　STM32VLDiscoveryに階乗、順列、組合せの計算をさせてみました。(STM32CubeIDE)
階乗は100!までを目標に。uint64_t(long long)で普通に計算すると直ぐに桁あふれするので、https://www.nak.ics.keio.ac.jp/class/hc/bignum.pdf を参考にさせて頂き、要素160個の配列で160桁までの計算を可能に。無限ループ中、その都度シリアル通信で数を文字列として与え、インタラクティブに計算させるようにもしました。細かいことですが桁が大きくなると見辛いので３桁ずつコンマ区切りで表示させる仕様にしました。

・階乗の結果は次のサイトと照合してみましたが正しく計算できているようです。http://www.suguru.jp/seminar/chuunyuusan/suukeisan/kaijou.html
101!も160桁に収まるので可能です。C言語では配列が渡せないので対処してあります。

void factorial(int k, int r, int a[keta]){
	int n, carry;
	a[0]=1;
	for(int i=1;i<keta;i++)a[i]=0;	//initialize a[]

	for(n=1+k-r;n<=k;n++){
		carry=0;
		for(int i=0;i<keta;i++){
			a[i]=a[i]*n;
			a[i]=a[i]+carry;
			carry=a[i]/10;
			a[i]=a[i]%10;
		}
	}

	int ketabegin=keta-1;
	n=ketabegin;
	while(1){
		if(a[n]!=0){ketabegin=n;break;}
		n--;
	}

	for(n=ketabegin;n>=0;n--){
		if((n+1)%3 ==0){
			if(n==ketabegin)printf("%d",a[n]);
			else printf(",%d",a[n]);
		}
		else printf("%d",a[n]);
	}
	printf("\r\n");
}

・順列は階乗の応用なので問題は特にありません。階乗や順列は配列の要素数次第です。

・組合せは、こちらで確認できます。https://keisan.casio.jp/exec/system/1161228812
再起呼び出しで行う 4.では34C17に5678秒ほど要するので34C17あたりが事実上限界です。3.では少し工夫して例えば
　　　10C6=10!/(6!x4!)
の計算を
　　　7x8x9x10/(1x2x3x4)=7/1x8/2x9/3x10/4
のようにしました。これでほぼ瞬時に計算してくれます。

/* 再帰呼び出し */
uint64_t combination(int n, int r){
	if (r == 0 || r == n) {return 1;
	} else if (r == 1) {   return n;}
	return combination(n - 1, r) + combination(n - 1, r - 1);
}

/* 組合せ改良版 */
uint64_t combination2(int n, int r){
	uint64_t c=1L;
	if(n-r>=r){
		for(int i=n-r+1;i<=n;i++){
			c *=i;
			c /=r-n+i;
		}
	}else{
		for(int i=r+1;i<=n;i++){
			c *=i;
			c /=i-r;
		}
	}
	return c;
}

/* 3桁ずつ整形して出力させる */
void ctocb( int cn, int cr, uint64_t val, int *a, int radix ){
	uint64_t v = val;//* 作業用(変換対象の値)
	int k;
	int ketae;
	int n = 0;//* 変換文字列の桁数記憶用
	while(v >= radix){//* 桁数を求める
		k=v % 10;
		a[n]=k;
		n++;
		v /= radix;
	}
	k=v % 10;
	a[n]=k;
	ketae=n+1;
	printf("%2dC%2d = ", cn, cr);
	for(int i=ketae;i>0;i--){
		if((i%3 ==1) && (i != 1)) printf("%d,",a[i-1]);
		else printf("%d",a[i-1]);
	}
}

問題はやはり桁あふれで63C32前後でおかしくなります。大きな数の割り算も配列を使ってできるそうなので、暇があれば改良したいと思います。
昔、PC9801DAフロッピーモデルでジーコジーコ計算させていたのが懐かしいです。
以下はMainです。

#define keta 160
int main(void){
	/* Reset of all peripherals, Initializes the Flash interface and the 
        Systick. */
	HAL_Init();
	/* Configure the system clock */
	SystemClock_Config();
	/* Initialize all configured peripherals */
	MX_GPIO_Init();
	MX_USART1_UART_Init();

	setbuf(stdout, NULL ); // printfを使用する前にバッファのリセット
	setbuf( stdin, NULL );
	printf("Math F, P, C Starts\r\n");
	printf("Help\r\n%s", HelpMsg1);
	printf("Input 1-5\r\n");
	uint64_t c=0L;
	int cb[54];
	int a[keta];
	char LINE[16];
	char* tmp;
	char* instr;
	char* str1;
	char* str2;
	int m,n,r;

	while (1){
		LINE[0] ='\0';         //* Flush Input Buffer */
		tmp='\0';
		instr='\0';
		str1='\0';
		str2='\0';
		m=0;
		n=0;
		r=0;

		printf(">" );
		if(!scanf( "%d",&m ))printf("\r\nIncorrect Command\r\n");
		if(m==5)break;	// END
		if(m>0 && m<5){
			switch(m){
			case 0:printf(">0 -> Help\r\n%s", HelpMsg1);
			case 1:
				setbuf(stdout, NULL );
				setbuf( stdin, NULL );
				printf(">%d. Factorial x(=<101 to 160 digits)\r\n",m);
				scanf( "%[^\r\n]",LINE );
				tmp = strdup(LINE);
				instr=trimString(tmp);
				if(!xatoi(&instr,&n))break;
				printf("%d!\r\n", n);
				factorial(n,n,a);
				break;
			case 2:
				setbuf(stdout, NULL );
				setbuf( stdin, NULL );
				printf(">%d. Permutation x,y : x(=<101)>y\r\n", m);
				scanf( "%[^\r\n]",LINE );
				tmp = strdup(LINE);
				instr=trimString(tmp);
				str1=strtok(instr, ",");
				str2=strtok(NULL,",");
				if(!xatoi(&str1,&n))break;
				if(!xatoi(&str2,&r))break;
				printf("%2dP%2d\r\n", n, r);
				factorial(n,r,a);
				break;
			case 3:
				setbuf(stdout, NULL );
				setbuf( stdin, NULL );
				printf(">%d. Combination2 x,y : x(=<62)>y\r\n", m);
				scanf( "%[^\r\n]",LINE );
				tmp = strdup(LINE);
				printf("\r\n");
				instr=trimString(tmp);
				str1=strtok(instr, ",");
				str2=strtok(NULL,",");
				if(!xatoi(&str1,&n))break;
				if(!xatoi(&str2,&r))break;
				c=combination2(n,r);
				ctocb(n,r,c,cb,10);
				printf("\r\n");
				break;
			case 4:
				setbuf(stdout, NULL );
				setbuf( stdin, NULL );
				printf(">%d. Combination x,y : x(=<33)>y\r\n", m);
				scanf( "%[^\r\n]",LINE );
				tmp = strdup(LINE);
				printf("\r\n");
				instr=trimString(tmp);
				str1=strtok(instr, ",");
				str2=strtok(NULL,",");
				if(!xatoi(&str1,&n))break;
				if(!xatoi(&str2,&r))break;
				Timer=0;
				c=combination(n,r);
				int took=Timer;
				ctocb(n,r,c,cb,10);
				printf(" ( %d sec)\r\n",took);
				break;
			default:break;
			}
		}
		printf("\r\n");
		fflush(stdin);	// 越えた文字列がstdin に残ってるので捨てる
	}
	printf(">5. END\r\n");
}