企业官网建设流程全解析

早前使用汇编语言写了一个高精度算法程式库, 意犹未尽, 今用C语言改写。C语言的好处是算法容易明白，修改也方便。此程式包含加、减、乘、除等基本算法，共有4个数组，分别叫a, b, c, d，各数组的长度由 SIZE 设定，最大可达 2730 字节，相当于10进制的 6570个小数位。

此程式共有5项输出。

第一项输出是 π, 采用马青公式（Machin Formula）公式：

先写出 calc_arctan 子程式，使用 divide_scalar 计算各分项比用 multiply 快很多。当计算时，因为除数= 239 x 239 = 57121 > 32767, 正好测试一下 divide_scalar 能否正确处理 remainder MSB 的进位。

// 任意精度 X*arctan(1/Y) 计算函数，结果保存在数组 a 中 // X, Y < 256 // 数组 b, c 被覆盖 void calc_arctan(unsigned int X, unsigned int Y) { unsigned int divisor = 3; unsigned char is_subtract = 1; // 1 = 减法, 0 = 加法（级数符号交替） // 计算第一项: term = X / Y clear(b); b[0] = X; divide_scalar(Y); // b = X / Y Y = Y*Y; // 下次除 Y^2 copy(b, a); // a = sum copy(b, c); // c = term // 迭代计算后续项 while (not_zero()) { copy(c, b); // b = term divide_scalar(Y); // b = term / Y copy(b, c); // c = 新 term divide_scalar(divisor); // b = 新 term / divisor if (is_subtract) { subtract(); // a = sum - 新 term } else { add(); // a = sum + 新 term } is_subtract = !is_subtract; // 切换加减符号 divisor += 2; } }

第二项输出是 ln 1.25，使用泰勒展开式 (Taylor series)，原公式收敛比较慢, 引入

把 1.25 转换为 1/9，由 multiply 改用 divide_scalar 计算各分项, 速度明显加快了。

即

​​​​​ 第三项输出是 e，公式比较简单,也是使用 divide_scalar，不断除 n, 然后n+1：

第四项输出是，使用牛顿迭代公式，每一次迭代 divide 都要执行 SIZE*8 次 shift_left 位元左移 (Bitwise Left Shift) 动作，速度最慢，计算1200位耗时 49 秒，计算6568位耗时30分鐘：

其中 ​​​​​。

唯一缩短运算时间的方法，就是把最耗时的 shift_left 子程式改为汇编语言。

// ASM 版本 bit shift_left() { #pragma ASM MOV DPTR, #SIZE-1 ; a 数组最后的地址 MOV WR4, #SIZE ; a 数组的大小 MOV A, R5 JZ CLR_C INC R4 CLR_C: CLR C LOOP: MOVX A, @DPTR RLC A ; MSB 移至 C MOVX @DPTR, A ; DEC DPTR DEC DPL MOV A, DPL XRL A, #0xFF ; XRL 不会影响 C JNZ LOOP1 DEC DPH LOOP1: DJNZ R5, LOOP ; DJNZ 不会影响 C DJNZ R4, LOOP #pragma ENDASM ; 返回 C }

使用ASM版本，divide 速度快了 30%。

最后一项输出是的平方，测试 multiply 的精确度，。

输出结果 (耗时1分14秒)：

Output pi up to 1200 decimal places: 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816096318595024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909216420198938095257201065485863278865936153381827968230301952035301852968995773622599413891249721775283479131515574857242454150695950829533116861727855889075098381754637464939319255060400927701671139009848824012 Output ln 1.25 up to 1200 decimal places: 0.223143551314209755766295090309834503374601085548007213671287872487391743768268333418407224100342235715963340980574191432352964757808415085568275114193553803690724495840440375272878778954558178115023454962871883866911484737848177562902024420172341248896755391915343869069834386977078762145959544090883810557672394536296476250152134418254417107481881140401651478372273686573452578249855026192763558094862502041388206129099704805174417608405611045487978507747746791146465972191457526488571334047924675817363189821622089684677155552892449421732245123818628048513040568976504515520642069147821498906237669983077674637864279162944840247509838362304317274098527186174441140506294263771766838297189401621233950728283407180817510452543533593598259516975775585290337904567271562482500085133416217789924759142508182505430778094286738537179002227654169816604980948324306350889362885190556380437265937652718605762558327948868241816974281891181560191589507098273620178149345023250074912708201734024432814516524393580452894575297017613473667720948338363167842157294002589126236390781103925354831318273259536107606622148483184797221222833426589156207184411570437132132079090449780503159557337609373173862093746778917 Output e up to 1200 decimal places: 2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427427466391932003059921817413596629043572900334295260595630738132328627943490763233829880753195251019011573834187930702154089149934884167509244761460668082264800168477411853742345442437107539077744992069551702761838606261331384583000752044933826560297606737113200709328709127443747047230696977209310141692836819025515108657463772111252389784425056953696770785449969967946864454905987931636889230098793127736178215424999229576351482208269895193668033182528869398496465105820939239829488793320362509443117301238197068416140397019837679320683282376464804295311802328782509819455815301756717361332069811250996181881593041690351598888519345807273866738589422879228499892086805825749279610484198444363463244968487560233624827041978623209002160990235304369941849146314093431738143640546253152096183690888707016768396424378140592714563549061303107208510383750510115747704171898610687396965521267154688957035035402123407849819334321068170121005627880235193033224745015853904730419957777093503660416997329725088687696640355570716226844716256079882651787134195124665201030592123667719432527867539855894489697096409 Output sqrt(2) up to 1200 decimal places: 1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641572735013846230912297024924836055850737212644121497099935831413222665927505592755799950501152782060571470109559971605970274534596862014728517418640889198609552329230484308714321450839762603627995251407989687253396546331808829640620615258352395054745750287759961729835575220337531857011354374603408498847160386899970699004815030544027790316454247823068492936918621580578463111596668713013015618568987237235288509264861249497715421833420428568606014682472077143585487415565706967765372022648544701585880162075847492265722600208558446652145839889394437092659180031138824646815708263010059485870400318648034219489727829064104507263688131373985525611732204024509122770022694112757362728049573810896750401836986836845072579936472906076299694138047565482372899718032680247442062926912485905218100445984215059112024944134172853147810580360337107730918286931471017111168391658172688941975871658215212822951848847208969463386289156288276595263514054226765323969461751129160240871551013515045538128756005263146801712740265396947024030051749531886292563138518816347800156936917688185237868405228783762938921430065586 Output 2 up to 1200 decimal places: 2.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

C 源码：​​

// 运行时 IRC 频率 24.000 MHz, 串口波特率 115200 bps #include <STC32G.h> #include <stdio.h> // 最大 SIZE 2730, DP 3284 #define SIZE 502 #define DP 600 // DP = (SIZE-1)*1.204, 输出 DP*2 个小数位 unsigned char xdata a[SIZE] = {0}; // 宣告 a 数组在 xdata unsigned char xdata b[SIZE] = {0}; // 宣告 b 数组在 xdata unsigned char xdata c[SIZE] = {0}; // 宣告 c 数组在 xdata unsigned char edata d[SIZE] = {0}; // 宣告 d 数组在 edata // 放在 code Flash RAM 内 unsigned char code hex2bcd[100] = { 0x00, 0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08, 0x09, 0x10, 0x11, 0x12, 0x13, 0x14, 0x15, 0x16, 0x17, 0x18, 0x19, 0x20, 0x21, 0x22, 0x23, 0x24, 0x25, 0x26, 0x27, 0x28, 0x29, 0x30, 0x31, 0x32, 0x33, 0x34, 0x35, 0x36, 0x37, 0x38, 0x39, 0x40, 0x41, 0x42, 0x43, 0x44, 0x45, 0x46, 0x47, 0x48, 0x49, 0x50, 0x51, 0x52, 0x53, 0x54, 0x55, 0x56, 0x57, 0x58, 0x59, 0x60, 0x61, 0x62, 0x63, 0x64, 0x65, 0x66, 0x67, 0x68, 0x69, 0x70, 0x71, 0x72, 0x73, 0x74, 0x75, 0x76, 0x77, 0x78, 0x79, 0x80, 0x81, 0x82, 0x83, 0x84, 0x85, 0x86, 0x87, 0x88, 0x89, 0x90, 0x91, 0x92, 0x93, 0x94, 0x95, 0x96, 0x97, 0x98, 0x99 }; void init_uart() { SCON = 0x50; // 串口1 模式1, 允许串口接收数据 T2L = 0xCC; // (65536 - (24000000UL / 115200) / 4) T2H = 0xFF; // (65536 - (24000000UL / 115200) / 4) >> 8 S1BRT = 1; // 选择定时器2 作为波特率发射器 T2x12 = 1; // 定时器 2 设置为 1T 模式（不分频） T2R = 1; // 启动定时器2 } // 透過串口1發送字符 void send(unsigned char dat) { SBUF = dat; while (!TI); // 等候完成發送 TI = 0; // Clear TI } // 透過串口1發送字串 void ptrstr(char *s) { while (*s) { send(*s); s++; } } // 发送两个 BCD 码至串口 void send_bcd(unsigned char x) { unsigned char ch; ch = (x >> 4) + '0'; send(ch); ch = (x & 0x0f) + '0'; send(ch); } // 将二进制数转为十进制 ASCII void send_bin2ascii(unsigned char val) { unsigned char h = val / 100; unsigned char t = (val / 10) % 10; unsigned char u = val % 10; if (h > 0) send((unsigned char) (h + '0')); if (h > 0 || t > 0) send((unsigned char) (t + '0')); send((unsigned char) (u + '0')); } // 任意精度乘法 a = a * q (只计算小数部分，保持 a[0] 整数部分不变) // 用于把 hex 转为 BCD unsigned int multiply_scalar(unsigned int q) { unsigned int i; unsigned long carry = 0; unsigned long xi; // 从最后一字节开始向前乘，直到索引 1（不包括索引 0） for (i = SIZE - 1; ; i--) { xi = (unsigned long) a[i] * q; xi += carry; carry = xi >> 8; a[i] = xi & 0xFF; if (i ==1) break; } return (unsigned int)carry; // 返回小数部分相乘后溢出到整数位的进位 } // 输出 a // a 数组被覆盖 void output() { unsigned int i; unsigned char dig; dig = a[0]; send_bin2ascii(dig); // 输出整数 send('.'); // 输出小数点 for (i = 0; i < DP; i++) { // 输出 DP*2 个小数位 dig = multiply_scalar(100); // a * 100转为10进制 send_bcd(hex2bcd[dig]); // 转为 BCD 发送 } } // 比较 a, b 大小 (若 a >= b 返回 1, 否则返回0) unsigned char is_greater() { unsigned int i; for (i = 0; i < SIZE; i++) { if (a[i] > b[i]) return 1; if (a[i] < b[i]) return 0; } return 1; // a=b } // 檢查 a 是否等于 b (a <> b 返回 1, 否则返回0) unsigned char diff() { unsigned int i; for (i = 0; i < SIZE; i++) { if (a[i] != b[i]) return 1; // 返回 1, x <> y } return 0; // 相等, 返回 0 } // 檢查 b = 0 unsigned char not_zero(){ unsigned int i; for (i = 0; i < SIZE; i++){ if (b[i] != 0) return 1; // 返回 1 = 不是零 } return 0; // 返回 0 = 零 } // a<< 1, a x 2, 向左移1位元 void shift_left(){ unsigned int i; // Global 1-bit shift left for (i = 0; i < SIZE - 1; i++) { a[i] = (a[i] << 1) | ((a[i+1] & 0x80) ? 1 : 0); } a[SIZE - 1] <<= 1; } // a >> 1, a / 2, 向右移1位元 void half(){ unsigned int i; // Global 1-bit shift right for (i = SIZE-1; i > 0; i--) { a[i] = (a[i] >> 1) | ((a[i-1] & 0x01) ? 0x80 : 0); } a[0] >>= 1; } // 清除数组 void clear(unsigned char x[]) { int i; for (i = 0; i < SIZE; i++) { x[i] = 0; } } // 复制 x 数组至 y 数组 void copy (unsigned char x[], unsigned char y[]) { unsigned int i; for (i = 0; i < SIZE; i++){ y[i] = x[i]; } } // a 数组与 b 数组交换 void swap(){ unsigned int i; unsigned char temp; for (i = 0; i < SIZE; i++){ temp = a[i]; a[i] = d[i]; d[i] = temp; } } // 任意精度加法: a = a + b void add() { unsigned int i; unsigned int carry = 0; unsigned int sum; for (i = SIZE - 1; ; i--) { sum = (unsigned int)a[i] + b[i] + carry; carry = sum >> 8; a[i] = sum & 0xFF; if (i == 0) break; } } // 任意精度減法: a = a - b void subtract() { unsigned int i; int borrow = 0; int sub; for (i = SIZE - 1; ; i--) { sub = (int)a[i] - b[i] - borrow; if (sub < 0) { sub += 256; borrow = 1; } else { borrow = 0; } a[i] = (unsigned char)sub; if (i == 0) break; } } // 任意精度乘法 c = a * b void multiply() { unsigned int i, j; unsigned long prod; // 清除 c 数组 (大小为 2 * SIZE) for (i = 0; i < (2 * SIZE); i++) { c[i] = 0; } // 带有快取的 256 进位小学教科书式乘法 for (i = SIZE - 1; ; i--) { unsigned char xi = a[i]; // 取被乘数的一个数位 a_i if (xi != 0) { // 跳过零 unsigned long carry = 0; for (j = SIZE - 1; ; j--) { // 依次去乘乘数的每一个数位 b_j unsigned long idx = (unsigned long)i + j; // 乘积的索引 prod = (unsigned long)xi * b[j] + c[idx] + carry; // 计算 积 = a_i * b_j + 进位 c[idx] = (unsigned char)(prod & 0xFF); // 存入乘积的低8位 carry = prod >> 8; // 进位是高8位 if (j == 0) { c[idx - 1] += (unsigned char)carry; // 最后的进位 break; } } } if (i == 0) break; } // --- 最后数位四捨五入 --- // 若最后数位 >= 0x80, 乘积的 LSB +1 if (c[SIZE] >= 0x80) { int carry_idx = SIZE - 1; while (carry_idx >= 0) { c[carry_idx]++; if (c[carry_idx] != 0) { break; // 沒有溢出, 停止 } carry_idx--; // 否则继续向前进1 } } } // 16位元任意精度除法: b = b / divisor void divide_scalar(unsigned int divisor) { unsigned int remainder = 0; int i=0; if (divisor == 0) return; // 跳过数组 b 前面的所有前导零字节 while (i < SIZE && b[i] == 0) { i++; } // 如果整个数组都是 0，则直接返回 if (i == SIZE) return; // 从第一个非零字节开始执行长除法 for (; i < SIZE; i++) { unsigned char quotient = 0; unsigned int bit_i; for (bit_i = 0; bit_i < 8; bit_i++) { unsigned char carry = (b[i] & 0x80) ? 1 : 0; unsigned char carry_in = (remainder & 0x8000) ? 1 : 0; b[i] <<= 1; // 若 divisor > 32767, remainder 有可能超过 65535, // 超过 unsigned int 范围, 进位会被舍弃, 需要用 carry_in 来保存进位 remainder = (remainder << 1) | carry; if (carry_in || (remainder >= divisor)) { remainder -= divisor; quotient |= (1 << (7 - bit_i)); } } b[i] = quotient; //更新商 } } // 任意精度除法: c = a / b // a 数组 被覆盖 void divide() { unsigned int byte_i; unsigned char bit_i; unsigned char quotient_byte = 0x00; // 整数部分, 不断相减 while (is_greater()) { // a > b subtract(); quotient_byte++; } c[0] = quotient_byte; // 分数部分 for (byte_i = 1; byte_i < SIZE; byte_i++) { quotient_byte = 0; for (bit_i = 0; bit_i < 8; bit_i++) { unsigned char carry_in = (a[0] & 0x80) ? 1 : 0; shift_left(); quotient_byte <<= 1; // 只有 carry_in = 0, 才进行比较 if (carry_in || is_greater()) { subtract(); quotient_byte |= 0x01; } } c[byte_i] = quotient_byte; } } // 任意精度 X*arctan(1/Y) 计算函数，结果保存在数组 a 中 // X, Y < 256 // b 数组, c 数组被覆盖 void calc_arctan(unsigned int X, unsigned int Y) { unsigned int divisor = 3; unsigned char is_subtract = 1; // 1 = 减法, 0 = 加法（级数符号交替） // 计算第一项: term = X / Y clear(b); b[0] = X; divide_scalar(Y); // b = X / Y Y = Y*Y; // 下次除 Y^2 copy(b, a); // a = sum copy(b, c); // c = term // 迭代计算后续项 while (not_zero()) { copy(c, b); // b = term divide_scalar(Y); // b = term / Y copy(b, c); // c = 新 term divide_scalar(divisor); // b = 新 term / divisor if (is_subtract) { subtract(); // a = sum - 新 term } else { add(); // a = sum + 新 term } is_subtract = !is_subtract; // 切换加减符号 divisor += 2; } } void main() { char buffer[40]; unsigned int divisor; WTST = 0; // 设置程序指令延时参数，赋值为0时速度设置为最快, 程序代码不等待 CKCON = 0; // 提高访问XRAM速度, xdata访问不增加等待 P_SW2 |= 0x80; // 扩展寄存器(XFR)访问使能 P3M0 = 0; // 串口1 输出/输入端口设置为准双向口 P3M1 = 0; init_uart(); // 运用马青公式求 pi: pi = 16 * arctan(1/5) - 4 * arctan(1/239) // 耗时 14 秒 (1200 位), 6分鐘 (6568 位) sprintf(buffer, "Output pi up to %d decimal places:\r\n", DP*2); ptrstr(buffer); //Calculate 4*arctan(1/239) calc_arctan(4, 239); swap(); // Calculate 16*arctan(1/5) calc_arctan(16, 5); copy(d, b); subtract(); output(); // 求 ln 1.25, (x-1)/(x+1) = 1/9 // 耗时 8 秒 (1200位), 3分36秒 (6568位) sprintf(buffer, "\r\nOutput ln 1.25 up to %d decimal places:\r\n", DP*2); ptrstr(buffer); clear(a); clear(b); b[0] = 1; divide_scalar(9); // b = term = 0.25/2.25 = 1/9 copy(b, a); // a = term shift_left(); // a = 2*term = sum copy(a, c); // c = 2*term = 新 term divisor = 3; while (not_zero()){ copy(c, b); // b = term divide_scalar(81); // b = term / 81 = 新 term copy(b, c); // c = 新 term divide_scalar(divisor); // b = 新 term / divisor add(); // a = a + b divisor+=2; } output(); // 求 e // 耗时 3 秒 (1200位), 1分26秒 (6568位) sprintf(buffer, "\r\nOutput e up to %d decimal places:\r\n", DP*2); ptrstr(buffer); clear(a); clear(b); a[0] = 2; // e = 2 b[0] = 1; divisor = 2; // 1/2 while (not_zero()){ divide_scalar(divisor); // b = 1 / n! add(); // a = a + 1/n! divisor++; // n+=1 } output(); // 求 sqrt(2) // 耗时 49 秒 (1200位), 30分鐘 (6568位) sprintf(buffer, "\r\nOutput sqrt(2) up to %d decimal places:\r\n", DP*2); ptrstr(buffer); clear(a); clear(b); a[0] = 1; while (diff()){ copy(a, b); // b = a(k-1) clear(a); a[0] =2; divide(); // c = 2/a(k-1) copy(c, a); add(); // a = 2/a(k-1) + a(k-1) half(); // a = (2/a(k-1) + a(k-1))/2 } copy(a, b); output(); sprintf(buffer, "\r\nOutput 2 up to %d decimal places:\r\n", DP*2); ptrstr(buffer); copy(b, a); multiply(); copy(c,a); output(); // a = 2 while(1); }