1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
|
program simple;
(* file exemples/pascal/simple.p: simple demo of Numerix
/*-----------------------------------------------------------------------+
| Copyright 2005-2006, Michel Quercia (michel.quercia@prepas.org) |
| |
| This file is part of Numerix. Numerix is free software; you can |
| redistribute it and/or modify it under the terms of the GNU Lesser |
| General Public License as published by the Free Software Foundation; |
| either version 2.1 of the License, or (at your option) any later |
| version. |
| |
| The Numerix Library is distributed in the hope that it will be |
| useful, but WITHOUT ANY WARRANTY; without even the implied warranty |
| of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
| Lesser General Public License for more details. |
| |
| You should have received a copy of the GNU Lesser General Public |
| License along with the GNU MP Library; see the file COPYING. If not, |
| write to the Free Software Foundation, Inc., 59 Temple Place - |
| Suite 330, Boston, MA 02111-1307, USA. |
+-----------------------------------------------------------------------+
| |
| Exemple d'utilisation de Numerix |
| |
+-----------------------------------------------------------------------*)
(* compute (sqrt(3) + sqrt(2))/(sqrt(3)-sqrt(2)) with n digits *)
uses _name_;
{$ifdef __GPC__}{$X+}{$endif}
(*
compute into res the quotient with n digits
res must be initialized by the caller
*)
procedure compute(var res:xint; n:int);
var a,b,d,d2,x,y:xint;
begin
(* d <- 10^n, d2 <- 10^(2n) *)
d := f_pow_1(5,n); shiftl(d,d,n);
d2 := f_sqr(d);
(* a <- round(sqrt(2*10^(2n+2))), b <- round(sqrt(3*10^(2n+2))) *)
a := f_mul_1(d2,200); gsqrt(a,a,1);
b := f_mul_1(d2,300); gsqrt(b,b,1);
(* res <- round(10^n*(b+a)/(b-a)) *)
x := f_add(b,a); mul(x,x,d);
y := f_sub(b,a);
gquo(res,x,y,1);
(* free temporary memory *)
xfree(d); xfree(d2);
xfree(a); xfree(b);
xfree(x); xfree(y);
end;
(* user interface *)
const test_n = 30;
const test_r = '9898979485566356196394568149411';
var i,n : int;
r : xint;
s : pchar;
test : boolean;
c : word;
begin
n := test_n;
r := xnew;
test:= false;
i := 1;
while i <= paramcount do begin
if paramstr(i) = '-test' then begin
test := true;
n := test_n;
i := paramcount+1;
end
else if (paramstr(i) = '-n') and (i < paramcount) then begin
val(paramstr(i+1),n,c);
c := c; (* to avoid a FPC warning about no-use of c *)
i := i+2;
end
else begin
writeln('usage: ',paramstr(0),' [-test] [-n <digits>]');
exit;
end
end;
compute(r,n);
s := string_of(r);
if test then begin
if strcmp(s,test_r) <> 0 then writeln('error in the ',paramstr(0),' test')
else writeln(paramstr(0),#9'test ok');
end
else writeln('r=',s);
(* done *)
strfree(s);
xfree(r);
end.
|
