1# This file is part of crlibm, the correctly rounded mathematical library,
2# which has been developed by the Arénaire project at École normale supérieure
3# de Lyon.
4#
5# Copyright (C) 2004-2011 David Defour, Catherine Daramy-Loirat,
6# Florent de Dinechin, Matthieu Gallet, Nicolas Gast, Christoph Quirin Lauter,
7# and Jean-Michel Muller
8#
9# This library is free software; you can redistribute it and/or
10# modify it under the terms of the GNU Lesser General Public
11# License as published by the Free Software Foundation; either
12# version 2.1 of the License, or (at your option) any later version.
13#
14# This library is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17# Lesser General Public License for more details.
18#
19# You should have received a copy of the GNU Lesser General Public
20# License along with this library; if not, write to the Free Software
21# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
22
23restart:
24Digits := 150;
25with (numapprox):with(orthopoly):
26interface(quiet=true);
27read "common-procedures.mpl";
28mkdir("TEMPCSH");
29
30
31#######################################################################
32# Some values to go to the .testdata
33# What is the first value that rounds to nearest to +inf ?
34ieeehexa(arcsinh(hexa2ieee(["7fefffff","ffffffff"])));
35
36# What is the first value that rounds to +inf to +inf ?
37
38
39######################################################################
40#First, some constants variables (used in the polynomial evaluations)
41n_double_ch := 11: # max degree for cosh's polynomial evaluation
42n_double_sh := 8:  # max degree for sinh's polynomial evaluation
43b_max := 2**(-9.): # max absolute input for polynomial evaluation
44
45######################################################################
46#some constants ...
47inv_ln_2 := 1/ln(2.):
48ln2_hi := hexa2ieee(["3FE62E42", "FEFA3800"]):
49ln2_lo := nearest(ln(2.)-ln2_hi):
50two_43_44 := 2^43 + 2^44:
51bias := convert(op(2,ieeehexa(ln(2.)/2.+two_43_44)),'decimal','hex'): #to get maximum index for the second range reduction ...
52
53
54######################################################################
55#Bounds of our evaluation
56max_input_ch := arccosh(hexa2ieee(["7FEFFFFF","FFFFFFFF"])):
57k_max_ch := ceil(max_input_ch / ln(2)):
58max_input_sh := arcsinh(hexa2ieee(["7FEFFFFF","FFFFFFFF"])):
59k_max_sh := ceil(max_input_sh / ln(2)):
60#towards +inf, we have cosh(x) = sinh(x) since exp(-x) is too small
61#so, we are sure of k_max_sh == k_max_ch
62k_max := k_max_ch;
63
64######################################################################
65# When can we ignore exp(-x) in front of exp(x) in the first step ?
66# We want the same error as in the general case
67k_max_csh_approx_exp := 35:
68tempxmax:=(k_max_csh_approx_exp-1)*log(2):
69eps_csh_approx_exp := exp(-tempxmax)/exp(tempxmax):
70log2(%);  #
71
72# When can we ignore exp(-x) in front of exp(x) in the second step ?
73# The worst case for exp for large arguments requires 115 bits
74k_max_csh_approx_exp_2 := 65:
75tempxmax:=(k_max_csh_approx_exp_2-1)*log(2):
76eps_csh_approx_exp_2 := exp(-tempxmax)/exp(tempxmax):
77log2(%); # 118
78
79
80######################################################################
81#The Taylor polynoms
82poly_ch :=series(cosh(x),x,n_double_ch):
83poly_ch := convert(poly_ch,polynom)-1;
84poly_sh :=series(sinh(x),x,n_double_sh):
85poly_sh := (convert(poly_sh,polynom))/x-1;
86
87
88
89
90
91
92####################################################################
93# secondary functions
94
95size_of_table := convert(op(2,ieeehexa(two_43_44+ln(2.)/2.)),decimal,hex);
96#returns the float which follow immediately the input
97next_float := proc(value)
98	local hex1,hex2,hexcat,result;
99	hex1:= op(1, value):
100	hex2:= op(2, value):
101	hexcat:= cat(hex1, hex2);
102	result := convert(convert(convert(hexcat,decimal,hex)+1+2**64,hex),string);
103	result := [substring(result,2..9), substring(result,10..18)];
104end:
105#compute the errors done in tabulated values for cosh
106delta_table_cosh_func := proc()
107	local result, i, value, temp, tmp, maxi;
108	value := ieeehexa(two_43_44-ln(2.)/2.);
109	result := 0;
110	maxi := 0;
111	for i from -size_of_table to size_of_table do
112		tmp := cosh(hexa2ieee(value)-two_43_44):
113		temp := nearest(tmp):
114		result := max(result, abs(tmp - temp - nearest(tmp-temp)));
115		maxi := max(maxi, abs(temp + nearest(tmp-temp)));
116		value:=next_float(value):
117	od:
118	result,maxi;
119end:
120#compute the error done in tabulated values for sinh
121delta_table_sinh_func := proc()
122	local result, i, value, temp, tmp, maxi;
123	value := ieeehexa(two_43_44-ln(2.)/2.);
124	result := 0;
125	maxi := 0;
126	for i from -size_of_table to size_of_table do
127		tmp := sinh(hexa2ieee(value)-two_43_44):
128		temp := nearest(tmp):
129		result := max(result, abs(tmp - temp - nearest(tmp-temp)));
130		maxi := max(maxi, abs(temp + nearest(tmp-temp)));
131		value:=next_float(value):
132	od:
133	result,maxi;
134end:
135
136
137#return the error on x_hi * y_lo
138Mul11_Error := proc(x,err_x, y, err_y)
139(2^(-53) * y + err_y) * (err_x + 1/2*ulp(x)) + x * err_y;
140end:
141
142
143#return the error on (x_hi * y_hi)_lo + x_lo * y_hi + x_hi * y_lo
144Mul43_Error := proc(x,err_x, y, err_y)
1451/2*ulp(3*2^(-53)* x * y) + 1/2*ulp(2*2^(-53) * x * y) + Mul11_Error(x, err_x, y, err_y) + Mul11_Error(y, err_y, x, err_x);
146end:
147
148
149
150cosh_0_35 := proc()
151local k, delta, maxi,epsilon, mini,cosh_0_35_max, delta_cosh_0_35:
152epsilon := 0; delta := 0:
153for k from -35 to -1 do
154delta_cosh_0_35 := 1/2*2^(-53)*ulp(1/(2^k)*(cosh_max + sinh_max)) + 1/(2^k)*(delta_sinh + delta_cosh):
155cosh_0_35_max := 1/(2^k)*(cosh_max + sinh_max):
156delta_cosh_0_35 :=  1/2*2^(-53)*ulp(cosh_0_35_max + 2^k * sinh_max) + delta_cosh_0_35 + 2^k  * delta_sinh:
157cosh_0_35_max := cosh_0_35_max + 2^k * sinh_max:
158delta_cosh_0_35 :=  1/2*2^(-53)*ulp(cosh_0_35_max + 2^k * cosh_max) + delta_cosh_0_35 + 2^k  * delta_cosh:
159cosh_0_35_max := cosh_0_35_max + 2^k * cosh_max:
160maxi := max(evalf(cosh((k-1/2)*ln(2))),evalf(cosh((k+1/2)*ln(2)))):
161mini := min(evalf(cosh((k-1/2)*ln(2))),evalf(cosh((k+1/2)*ln(2)))):
162
163epsilon := max(epsilon, delta_cosh_0_35/mini):
164delta := max(delta, delta_cosh_0_35):
165od;
166for k from 1 to 35 do
167delta_cosh_0_35 := 1/2*2^(-53)*ulp(1/(2^k)*(cosh_max + sinh_max)) + 1/(2^k)*(delta_sinh + delta_cosh):
168cosh_0_35_max := 1/(2^k)*(cosh_max + sinh_max):
169delta_cosh_0_35 :=  1/2*2^(-53)*ulp(cosh_0_35_max + 2^k * sinh_max) + delta_cosh_0_35 + 2^k  * delta_sinh:
170cosh_0_35_max := cosh_0_35_max + 2^k * sinh_max:
171delta_cosh_0_35 :=  1/2*2^(-53)*ulp(cosh_0_35_max + 2^k * cosh_max) + delta_cosh_0_35 + 2^k  * delta_cosh:
172cosh_0_35_max := cosh_0_35_max + 2^k * cosh_max:
173
174maxi := max(evalf(cosh((k-1/2)*ln(2))),evalf(cosh((k+1/2)*ln(2)))):
175mini := min(evalf(cosh((k-1/2)*ln(2))),evalf(cosh((k+1/2)*ln(2)))):
176epsilon := max(epsilon, delta_cosh_0_35/mini):
177delta := max(delta, delta_cosh_0_35):
178od;
179delta, epsilon;
180end:
181
182
183cosh_35_inf := proc()
184local k, delta, maxi,epsilon, mini,cosh_0_35_max, delta_cosh_0_35:
185epsilon := 0; delta := 0:
186for k from 35 to 1025 do
187delta_cosh_0_35 := 1/2*2^(-53)*ulp(2^k*(cosh_max + sinh_max)) + 2^k*(delta_sinh + delta_cosh) + 1/(2^k)*(cosh_max + sinh_max + delta_sinh + delta_cosh):
188cosh_0_35_max := 2^k*(cosh_max + sinh_max):
189
190maxi := max(evalf(cosh((k-1/2)*ln(2))),evalf(cosh((k+1/2)*ln(2)))):
191mini := min(evalf(cosh((k-1/2)*ln(2))),evalf(cosh((k+1/2)*ln(2)))):
192epsilon := max(epsilon, delta_cosh_0_35/mini):
193delta := max(delta, delta_cosh_0_35):
194od;
195delta, epsilon;
196end:
197#############################################################
198
199
200
201
202
203
204
205
206######################################################################
207#now we can begin the proof
208
209######################################################################
210#first, we must compute the error created by the first range reduction
211# CODY and WAITE  Argument reduction
212
213ln2 := ln(2.);
214invln2:= nearest(1/ln2);
215reminvln2 := evalf(1/ln2 - invln2);
216expln2:=ieeedouble(ln2)[2]:   #get the exponent of ln2 in its IEEE representation
217
218bits_ln2_hi_0 := ceil(log2(k_max));
219# 1/2 <= ln2/2^(expln2 + 1) <  1
220ln2_hi := round(evalf(ln2 * 2^(52 - bits_ln2_hi_0 - expln2)))/2^(52 - bits_ln2_hi_0 - expln2);#this trick is used to get a truncated mantissa for ln2_hi
221#ln2_hi is a now exactly representable in the IEEE format and bits_ln2_hi_0 last bits of its mantissa are set to 0
222#and bits_ln2_hi_0 is exactly the max number of bits to represent k
223#so the k * ln2_hi-product is exact :)
224ln2_lo:=nearest(ln2 - ln2_hi):
225
226# The error in this case (we need absolute error)
227delta_repr_ln2 := abs(ln2 - ln2_hi - ln2_lo);
228delta_round := evalf(1/2 * ulp(ln2_lo));
229delta_cody_waite := k_max * (delta_repr_ln2 + delta_round);
230
231delta_b := delta_repr_ln2 + delta_round + delta_cody_waite;
232
233#we have 2 cases:
234#  * k != 0 and we have an inexact range reduction
235#  * k == 0 and we have no range reduction and then delta_range_reduc = 0
236
237#the second range reduction is exact, so it doesn't introduce new error
238#after this second range reduction, we have a argument <= 2^(-9)
239#'mathematical' reductions:
240# x = k * ln(2) + y
241# y = a + b
242#'true' reductions:
243# x = k * (ln2_hi + ln2_lo) + (b_hi + b_lo) + table_index_float
244#   with table_index_float = table_index * 2^(-8)
245
246#we'll use the following mathematical formulaes :
247#cosh(a + b) = cosh(a) * cosh(b) + sinh(a) * sinh(b)
248#sinh(a + b) = sinh(a) * cosh(b) + sinh(b) * cosh(a)
249#sinh(a) and cosh(a) are tabulated as double double
250# and we use Taylor series to compute approximations of the sums
251
252#computation of the absolute error in the tabulated values :
253delta_ca := delta_table_cosh_func()[1];
254delta_sa := delta_table_sinh_func()[1];
255ca_max := delta_table_cosh_func()[2];
256sa_max := delta_table_sinh_func()[2];
257
258
259#now we must compute the error done in polynomial evaluation
260#we use   cosh(b) = 1 + sum(b^(2*k)/(2*k!), k > 0)
261#         sinh(b) = b * (1 + sum(b^(2*k)/(2*k+1!), k > 0))
262
263#both used polynoms are even (x�, x^4, x^6, ...) and we can use this fact
264y_max := b_max ^ 2;
265delta_y := 1/2*ulp(y_max) + delta_b^2;
266# remove the first x and compute the polynomial of y = x�
267poly_ch2 :=  subs(x=sqrt(y), expand(poly_ch));
268poly_sh2 :=  subs(x=sqrt(y), expand(poly_sh));
269errlist_cosh := errlist_quickphase_horner(degree(poly_ch2), 0, 0, 2^(-53), 2^(-70)):
270errlist_sinh := errlist_quickphase_horner(degree(poly_sh2), 0, 0, 2^(-53), 2^(-70)):
271
272#error between effective result and theorical polynomial result
273rounding_error_tcb := compute_horner_rounding_error(poly_ch2, y, y_max, errlist_cosh, true);
274rounding_error_tsb := compute_horner_rounding_error(poly_sh2, y, y_max, errlist_sinh, true);
275
276#error between therical polynomial result and cosh value
277approx_error_tcb := infnorm((cosh(x)-1-poly_ch)/(cosh(x)-1),x= -b_max..b_max);
278approx_error_tsb := infnorm((sinh(x)/x-1-poly_sh)/(sinh(x)/x-1),x= -b_max..b_max);
279
280delta_tcb := rounding_error_tcb[2] + approx_error_tcb;
281delta_tsb := rounding_error_tsb[2] + approx_error_tsb;
282tsb_max := rounding_error_tsb[4];
283tcb_max := rounding_error_tcb[4];
284
285#now we must do the first reconstruction, which correspond to cosh(a + b) = cosh(a) * cosh(b) + sinh(a) * sinh(b)
286#first case : sinh(a) = 0 = sa (= sa_hi + sa_lo in the C code)
287#             cosh(a) = 1 = ca (= ca_hi + ca_lo in the C code)
288#there is no error on sa and caj
289delta_cosh0 := delta_tcb;
290cosh0_max := 1+tcb_max;
291delta_sinh0 := delta_b:
292sinh0_max := b_max:
293delta_sinh0 := delta_sinh0 + (tsb_max + delta_tsb)*(b_max + delta_b + 1/2*ulp(b_max)) - b_max * tsb_max + 1/2*ulp(sinh0_max + tsb_max*b_max):
294sinh0_max := sinh0_max + tsb_max*b_max:
295delta_sinh0 := delta_sinh0 + 2^(-53)*1/2*ulp(sinh0_max + tsb_max*b_max);
296sinh0_max := sinh0_max + tsb_max*b_max;
297
298#second case : sinh(a) <> 0
299#there is a delta_table_cosh and delta_table_sinh absolute error on ca and sa.
300
301delta_cosh1 := delta_ca:
302cosh1_max := 2^(-53)*ca_max:
303delta_cosh1 := 1/2*ulp(cosh1_max + 3*2^(-53)*b_max * sa_max) + delta_cosh1 + Mul43_Error(b_max, delta_b, sa_max, delta_sa):
304cosh1_max := cosh1_max + 3 * 2^(-53) * b_max * sa_max:
305delta_cosh1 := 1/2*ulp(cosh1_max + b_max * sa_max * tsb_max) + delta_cosh1 + ((sa_max+1/2*ulp(sa_max)+delta_sa)*(b_max+ulp(b_max)+delta_b)*(tsb_max+delta_tsb)-tsb_max*sa_max*b_max):
306cosh1_max := cosh1_max + b_max * sa_max * tsb_max:
307delta_cosh1 := 1/2*ulp(cosh1_max + tcb_max*ca_max) + delta_cosh1 + ((ca_max + 1/2*ulp(ca_max) + delta_ca)*(tcb_max-delta_tcb)-ca_max*tcb_max):
308cosh11_max := cosh1_max + tcb_max*ca_max:
309delta_cosh1 := 1/2*ulp(cosh1_max + b_max * sa_max) + delta_cosh1:
310cosh11_max := cosh1_max + b_max * sa_max:
311delta_cosh1 := 1/2*2^(-53)*ulp(cosh1_max + ca_max) + delta_cosh1;
312cosh11_max := cosh1_max + ca_max;
313
314cosh_max := max(cosh0_max, cosh1_max);
315delta_cosh := max(delta_cosh0, delta_cosh1);
316
317delta_sinh1 := delta_sa:
318sinh1_max := 2^(-53)*sa_max:
319delta_sinh1 := 1/2*ulp(sinh1_max + 3*2^(-53)*ca_max * b_max) + delta_sinh1 + Mul43_Error(ca_max, delta_ca, b_max, delta_b):
320sinh1_max := sinh1_max + 3*2^(-53)*ca_max * b_max:
321delta_sinh1 := delta_sinh1 + 1/2*ulp(sinh1_max + sa_max * tcb_max):
322sinh1_max := sinh1_max + sa_max * tcb_max:
323delta_sinh1 := 1/2*ulp(sinh1_max + b_max * ca_max * tsb_max) + delta_sinh1 + ((ca_max+1/2*ulp(ca_max)+delta_ca)*(b_max+ulp(b_max)+delta_b)*(tsb_max+delta_tsb)-tsb_max*ca_max*b_max):
324sinh1_max := sinh1_max + b_max * sa_max * tsb_max:
325delta_sinh1 := delta_sinh1 + 1/2*2^(-53)*ulp(sinh1_max + b_max * ca_max):
326sinh1_max := sinh1_max + b_max * ca_max:
327delta_sinh1 := delta_sinh1 + 2^(-53)*1/2*ulp(sinh1_max + sa_max);
328sinh1_max := sinh1_max + sa_max;
329
330sinh_max := max(sinh1_max, sinh0_max);
331delta_sinh := max(delta_sinh0, delta_sinh1);
332#so we have the error done on cosh(y) and sinh(y)
333#now we must compute the error done on the last reconstruction
334#there are many cases
335#we begin by 0 < |k| < 35
336epsilon_cosh_0_35 := cosh_0_35()[2];
337#|k| > 35
338epsilon_cosh_35_inf := cosh_35_inf()[2];
339
340#rounding constant
341maxepsilon_csh := max(epsilon_cosh_35_inf, epsilon_cosh_0_35, delta_cosh):
342round_cst_csh := evalf(compute_rn_constant(maxepsilon_csh));;
343
344
345
346#################################################################################################"
347#now some functions used to build the .h header file.
348#################################################################################################"
349IEEE2db_number_BE := proc(ieee_number)
350local hex1, hex2, hexcat;
351hexcat=ieeehexa(ieee_number);
352hex1:= op(1, ieeehexa(ieee_number)):
353hex2:= op(2, ieeehexa(ieee_number)):
354cat(cat("{{0x"||hex1||",0x"||hex2||"}};     /*",sprintf("%.10e",ieee_number)),"*/  \n");
355end:
356IEEE2db_number_LE := proc(ieee_number)
357local hex1, hex2, hexcat;
358hexcat=ieeehexa(ieee_number);
359hex1:= op(2, ieeehexa(ieee_number)):
360hex2:= op(1, ieeehexa(ieee_number)):
361cat(cat("{{0x"||hex1||",0x"||hex2||"}};     /*",sprintf("%.10e",ieee_number)),"*/  \n");
362end:
363IEEE2db_number := proc(ieee_number, big_little)
364if (big_little = 1) then
365IEEE2db_number_BE(ieee_number):
366else
367IEEE2db_number_LE(ieee_number):
368fi;
369end:
370lo_part := proc(ieee_number)
371if (ieee_number <> 0) then
372ieee_number - nearest(ieee_number);
373else
3740;
375end if;
376end:
377IEEE2db_db_number_sinh := proc(number,big_little)
378local hexstring1,hexstring2, hex1,hex2,hex3,hex4;
379hexstring1 := ieeehexa(number);
380hexstring2 := ieeehexa(lo_part(number));
381if (big_little = 1) then
382hex1:= op(1, hexstring1):
383hex2:= op(2, hexstring1):
384hex3:= op(1, hexstring2):
385hex4:= op(2, hexstring2):
386else
387hex1:= op(2, hexstring1):
388hex2:= op(1, hexstring1):
389hex3:= op(2, hexstring2):
390hex4:= op(1, hexstring2):
391fi;
392"	  {{0x"||hex1||", 0x"||hex2||"}}, {{0x"||hex3||", 0x"||hex4||"}}},\n";
393end:
394
395IEEE2db_db_number_cosh := proc(number,big_little)
396local hexstring1,hexstring2, hex1,hex2,hex3,hex4;
397hexstring1 := ieeehexa(number);
398hexstring2 := ieeehexa(lo_part(number));
399if (big_little = 1) then
400hex1:= op(1, hexstring1):
401hex2:= op(2, hexstring1):
402hex3:= op(1, hexstring2):
403hex4:= op(2, hexstring2):
404else
405hex1:= op(2, hexstring1):
406hex2:= op(1, hexstring1):
407hex3:= op(2, hexstring2):
408hex4:= op(1, hexstring2):
409fi;
410"	{{{0x"||hex1||", 0x"||hex2||"}}, {{0x"||hex3||", 0x"||hex4||"}},\n";
411end:
412
413#####################################################################################
414
415
416
417#now, we can produce the header file !
418round_cst_cosh := 1.0020:
419round_cst_sinh := round_cst_cosh:
420filename := "TEMPCSH/csh_fast.h":
421fd := fopen(filename, WRITE, TEXT):
422fprintf(fd, "\n /* File generated by maple/csh.mpl */ \n"):
423fprintf(fd, "\n"):
424fprintf(fd, " static double maxepsilon_csh = %1.30e ;\n", maxepsilon_csh):
425fprintf(fd, " static double round_cst_csh  = %1.30e ;\n", round_cst_csh):
426fprintf(fd, "\n"):
427
428
429fprintf(fd, "#ifdef WORDS_BIGENDIAN  \n"):
430
431for big_little from 1 to 2 do
432	if (big_little = 2) then
433		fprintf(fd,                "#else  \n"):
434	fi:
435	fprintf(fd, cat(  "  static db_number const inv_ln_2 =     ",IEEE2db_number(inv_ln_2,big_little))):
436	fprintf(fd, cat(  "  static db_number const ln2_hi =       ",IEEE2db_number(ln2_hi,big_little))):
437	fprintf(fd, cat(  "  static db_number const ln2_lo =       ",IEEE2db_number(ln2_lo,big_little))):
438	fprintf(fd, cat(  "  static db_number const two_43_44 =    ", IEEE2db_number(two_43_44,big_little))):
439	fprintf(fd, cat(  "  static db_number const two_minus_30 = ", IEEE2db_number(2**(-40),big_little))):
440	fprintf(fd,       "  static int const bias = %d ;\n", bias):
441	fprintf(fd, "\n");
442
443	fprintf(fd,"/* some bounds */ \n"):
444	fprintf(fd, cat(  "  static db_number const max_input_csh =  ",IEEE2db_number(max_input_ch,big_little))):
445
446    fprintf(fd, "\n"):
447	fprintf(fd, cat(cat("  static const db_number cosh_sinh_table[",convert(2*size_of_table+1,string)),"][4] = { \n"));
448	vvalue := ieeehexa(two_43_44-ln(2.)/2.);
449	for i from -size_of_table to size_of_table do
450		fprintf(fd, IEEE2db_db_number_cosh(cosh(hexa2ieee(vvalue)-two_43_44),big_little));
451		fprintf(fd, IEEE2db_db_number_sinh(sinh(hexa2ieee(vvalue)-two_43_44),big_little));
452		vvalue:=next_float(vvalue):
453	od:
454	fprintf(fd,"}; \n");
455
456	fprintf(fd,"/* the coefficients for the cosh-approximations */ \n"):
457	for i from 1 to (n_double_ch/2) do
458		fprintf(fd, cat(cat(cat(  "  static const db_number c",convert(2*(i-1),string))," =   "),IEEE2db_number(coeff(poly_ch+1,x,2*(i-1)),big_little))):
459	od:
460
461	fprintf(fd,"/* the coefficients for the sinh-approximations */\n"):
462	for i from 1 to (n_double_sh/2) do
463		fprintf(fd, cat(cat(cat(  "  static  const db_number s",convert(2*i-1,string))," =   "),IEEE2db_number(coeff(poly_sh+1,x,2*(i-1)),big_little))):
464	od:
465od:
466fprintf(fd, "#endif  \n"):
467
468fclose(fd):
469