ubuntu-fr Scripts - Math

Merge lp://staging/~ubuntu-fr-scripts/ufrs-math/lib-bc into lp://staging/ufrs-math

  1. lib-bc
  2. Merge into math
Proposed by Jean-Mi
Status: Merged
Approved by: Jean-Mi
Approved revision: 4
Merge reported by: Jean-Mi
Merged at revision: not available
Proposed branch: lp://staging/~ubuntu-fr-scripts/ufrs-math/lib-bc
Merge into: lp://staging/ufrs-math
Diff against target: None lines
To merge this branch: bzr merge lp://staging/~ubuntu-fr-scripts/ufrs-math/lib-bc
Reviewer Review Type Date Requested Status
Jean-Mi Approve
Review via email: mp+5058@code.staging.launchpad.net
To post a comment you must log in.
Revision history for this message
Jean-Mi (jeanmi) wrote :

Le répertoire est prêt pour l'intégration dans la branche principale

Revision history for this message
Jean-Mi (jeanmi) wrote :

ok c'est bon

review: Approve

Preview Diff

[H/L] Next/Prev Comment, [J/K] Next/Prev File, [N/P] Next/Prev Hunk
1=== added directory 'bc-lib'
2=== added file 'bc-lib/README.bc-lib'
3--- bc-lib/README.bc-lib 1970-01-01 00:00:00 +0000
4+++ bc-lib/README.bc-lib 2009-03-31 00:25:36 +0000
5@@ -0,0 +1,10 @@
6+# Bibliothèque de fonctions pour GNU bc
7+#
8+# Auteurs : Jean-Michel Juzan ainsi que d'autre membres
9+# de la communauté du libre
10+#
11+# export BCLIB="/usr/share/ufrs-math/bc-lib"
12+# Usage : bc -l $BCLIB/nom_bibliothèque.bc
13+#
14+# Dépend de bc
15+#
16
17=== added file 'bc-lib/anglelog.bc'
18--- bc-lib/anglelog.bc 1970-01-01 00:00:00 +0000
19+++ bc-lib/anglelog.bc 2009-03-21 02:39:06 +0000
20@@ -0,0 +1,36 @@
21+### AngleLog.BC - logarithmic function that are semi-linear,
22+# changing gradient at each power of the base for GNU bc.
23+
24+# run with funcs.bc
25+
26+## Specific functions for base 10
27+
28+# Generate the x-th number in the sequence:
29+# 1,2,3,4,5,6,7,8,9,10,20,30,...,90,100,200,300,...,900,1000,2000,...
30+# e.g. al10(0) = 1 ; al10(11) = 30
31+define al10(x) {
32+ return( (rem(x,9)+1) * 10^int(x/9) )
33+}
34+
35+# Invert the above
36+# e.g. ial10(50) = 13 ; ial10(1) = 0
37+define ial10(x) {
38+ auto k;
39+ k = int(log(10,x))
40+ return( 9*k + x/10^k - 1)
41+}
42+
43+## Generalised functions for any base b
44+
45+# al(10,x) is equivalent to the above al10(x)
46+define al(b,x) {
47+ return( (rem(x,b-1)+1) * b^int(x/(b-1)) )
48+}
49+
50+# Invert the previous function
51+define ial(b,x) {
52+ auto k
53+ k = int(log(b,x))
54+ return( (b-1)*k + x/b^k - 1)
55+}
56+
57
58=== added file 'bc-lib/cf.bc'
59--- bc-lib/cf.bc 1970-01-01 00:00:00 +0000
60+++ bc-lib/cf.bc 2009-03-21 02:39:06 +0000
61@@ -0,0 +1,97 @@
62+### CF.BC - Continued fraction experimentation for GNU bc
63+
64+# Create a continued fraction representation of x in the global array cf[]
65+define cf(x) {
66+ auto sign, i, s3
67+ s3=int(scale/2);
68+ sign=sgn(x);x*=sign
69+ for(i=0;i<=s3;i++) {
70+ cf[i]=int(x)
71+ if(cf[i]>10^s3)break;
72+ x-=cf[i]
73+ if(x==0){i+=1;break}
74+ x=1/x
75+ }
76+ cf[i]=0
77+}
78+
79+# Create a continued fraction representation of x in the global array cf[]
80+# using alternating signs
81+define cf2(x) {
82+ auto sign, i, s3
83+ s3=int(scale/2);
84+ sign=sgn(x);x*=sign
85+ for(i=0;i<=s3;i++) {
86+ cf[i]=1+int(x)
87+ if(cf[i]>10^s3)break;
88+ x-=cf[i]
89+ if(i!=2*int(i/2))cf[i]*=-1
90+ if(x==0){i+=1;break}
91+ x=1/-x
92+ }
93+ cf[i]=0
94+}
95+
96+# Output the global array cf[] formatted as a continued fraction
97+define printcf() {
98+ auto i;
99+ print "[", cf[0], "; ";
100+ i=1;if(cf[i]!=0)for(i=1;cf[i+1]!=0;i++)print cf[i], ", ";
101+ print cf[i], "]\n";
102+}
103+
104+# Convert global array cf[] back into a number
105+define cf2num() {
106+ auto n, i;
107+ for(i=0;cf[i]!=0;i++){}
108+ n=cf[--i]
109+ for(;i>=0;i--)n=1/n+cf[i]
110+ return(n);
111+}
112+
113+# Turn the alternating signed continued fraction representation of x back
114+# into a number by first taking the absolute value of the chain
115+define cfflip(x) {
116+ cf2(x)
117+ for(i=0;cf[i]!=0;i++)cf[i]=abs(cf[i])
118+ return(cf2num());
119+}
120+
121+# Turn the alternating signed continued fraction representation of x back
122+# into a number by first taking the absolute value less one of the chain
123+define cfflip1(x) {
124+ cf2(x)
125+ for(i=0;cf[i]!=0;i++)cf[i]=abs(cf[i])-1
126+ return(cf2num());
127+}
128+
129+# Turn the binary representation of x into a continued fraction-like
130+# structure using global array cf[].
131+# e.g. 0.1001000011111101110111 -> [0;1,2,1,4,6,1,3,1,3]
132+define bincf(x) {
133+ auto i,b,bb,n,j;
134+ x=abs(x);if(x>1)x=frac(x)
135+ x*=2;b=int(x);x-=b;n=1;j=1;cf[0]=0
136+ for(i=0;i<scale;i++){
137+ x*=2;bb=int(x)
138+ if(bb==b){n+=1}else{cf[j]=n;j+=1;n=1}
139+ b=bb;x-=b
140+ }
141+ if(n!=0){cf[j]=n;j+=1}
142+ cf[j-1]*=b
143+}
144+
145+# Read a continued fraction as a representation of alternating groups of
146+# bits in a binary number.
147+# e.g. [1;4,1,4,1,4,1,4,1...] -> 0.11110111101111011110...
148+define cfbin() {
149+ auto n,i,b,x;
150+ #cf[0]=0;
151+ for(i=1;cf[i]!=0;i++){}
152+ n=cf[--i];b=1;x=0
153+ for(;i>=0;i--){
154+ for(j=0;j<n;j++)x=x/2+b
155+ b=1-b
156+ }
157+ return(x)
158+}
159
160=== added file 'bc-lib/collatz.bc'
161--- bc-lib/collatz.bc 1970-01-01 00:00:00 +0000
162+++ bc-lib/collatz.bc 2009-03-21 02:39:06 +0000
163@@ -0,0 +1,113 @@
164+### Collatz.BC - The 3x+1 or hailstones problem for GNU bc
165+
166+## Here we'll be using the modified rules of
167+## Odd -> (3x+1)/2
168+## Even -> x/2
169+## ...as otherwise the usual odd step (3x+1) is always followed by
170+## the even step. This way each rule has a 50/50 distribtion.
171+
172+# Generate the next hailstone 5 -> (3*5+1)/2 = 8; 8 -> 8/2 = 4;
173+define cz_next(x) {
174+ auto os;if(x<=1)return(1)
175+ os=scale;scale=0;x/=1
176+ if(x%2)x=3*x+1;x/=2
177+ scale=os
178+ return(x)
179+}
180+
181+# Take a guess at the previous hailstone - since in some cases there are
182+# two choices, this function always chooses the lower option.
183+# e.g. 5 comes from both 3 [(3*(3)+1)/2] and 10 [10/2]
184+# - this function chooses 3.
185+define cz_prev(x) {
186+ auto os;if(x<=1)return(1)
187+ os=scale;scale=0;x/=1
188+ x=2*x-1;if(x%3){x+=1}else{x/=3}
189+ scale=os
190+ return(x)
191+}
192+
193+# Determine the hailstone chain length to 1
194+# e.g. [5 -> 8 -> 4 -> 2 -> 1] = 4 steps
195+define cz_chlen(x) {
196+ auto os,c;if(x<=1)return(0)
197+ os=scale;scale=0;x/=1
198+ c=0
199+ while(x!=1){
200+ c+=1
201+ if(x%2)x=3*x+1;x/=2
202+ }
203+ scale=os
204+ return(c)
205+}
206+
207+# Add the values of all the hailstones in a chain.
208+# [5 -> 8 -> 4 -> 2 -> 1] => 5 + 8 + 4 + 2 + 1 = 20
209+define cz_chsum(x) {
210+ auto os,s;if(x<=1)return(0)
211+ os=scale;scale=0;x/=1
212+ s=x
213+ while(x!=1){
214+ if(x%2)x=3*x+1;x/=2
215+ s+=x
216+ }
217+ scale=os
218+ return(s)
219+}
220+
221+# Print the whole hailstone chain of numbers from x
222+define cz_chain(x) {
223+ auto os;if(x<=1)return(0)
224+ os=scale;scale=0;x/=1
225+ while(x!=1){
226+ x
227+ if(x%2)x=3*x+1;x/=2
228+ }
229+ scale=os
230+ return(1)
231+}
232+
233+# Find the highest 'altitude' that the hailstone chain for x reaches
234+define cz_chmax(x) {
235+ auto os,m;if(x<=1)return(0)
236+ os=scale;scale=0;x/=1
237+ m=0
238+ while(x!=1){
239+ if(x>m)m=x
240+ if(x%2)x=3*x+1;x/=2
241+ }
242+ scale=os
243+ return(m)
244+}
245+
246+# Generate a binary number representing a hailstone choice path
247+define cz_8tp(x) { # 8-tree path
248+# all numbers lead back to 8, this function maps the path as a binary
249+# number with the bit representing 8 as most significant.
250+# (8) = 1, (5 -> 8) = 10, (16 -> 8) = 11, (10 -> 5 -> 8) = 101, etc
251+ auto os,p,t;if(x<4)return(0)
252+ os=scale;scale=0;x/=1
253+ p=0;t=1
254+ while(x!=8){
255+ if(x%2){x=(3*x+1)/2}else{x/=2;p+=t}
256+ t*=2
257+ }
258+ p+=t
259+ scale=os
260+ return(p)
261+}
262+
263+# Generate a hailstone from its path's binary number
264+define cz_i8tp(p) { # convert an 8-tree path back into a number
265+ auto os,msb,x;if(p<1)return(0)
266+ os=scale;scale=0;p/=1
267+ msb=1;while(msb<=p)msb*=2;msb/=2
268+ x=8;p-=msb;msb/=2
269+ while(msb){
270+ while(msb>p){scale=os;x=(2*x-1)/3;scale=0; msb/=2}
271+ if(msb){scale=os;x*=2;scale=0; p-=msb;msb/=2}
272+ }
273+ p=x/1;if(p==x)x=p # remove zeroes after dp if there are any
274+ scale=os
275+ return(x)
276+}
277
278=== added file 'bc-lib/complex.bc'
279--- bc-lib/complex.bc 1970-01-01 00:00:00 +0000
280+++ bc-lib/complex.bc 2009-03-21 02:39:06 +0000
281@@ -0,0 +1,200 @@
282+### Complex.BC - Rudimentary complex number handling for GNU BC
283+
284+## Not to be regarded as suitable for any purpose
285+## Not guaranteed to return correct answers
286+
287+# Read a digit from the decimal representation of a number
288+# read_digit(3210.987, -1) returns 9 as that is the digit in
289+# the 10^-1 position.
290+define read_digit(number, place) {
291+ auto os,d;os = scale
292+ number*=10^-place
293+ scale=0
294+ d = (number%10)/1
295+ scale=os
296+ return(d)
297+}
298+
299+# MakeComplex: Turn two numbers into one by interleaving the digits
300+define makecomplex(r,i) {
301+ auto c, rd, id, ri, ii, j, sign
302+
303+ r /= 1
304+ i /= 1
305+ ri = length(r)-scale
306+ ii = length(i)-scale
307+
308+ sign = 1
309+ if (r<0){r*=-1;sign+=2}
310+ if (i<0){i*=-1;sign+=1}
311+
312+ if(ri<ii)ri=ii
313+ scale *= 2
314+
315+ c = 0/1
316+
317+ for(j=-scale;j<ri;j++) {
318+ rd = read_digit(r, j)
319+ id = read_digit(i, j)
320+ c += (rd*10+id)*100^j
321+ }
322+
323+ c += 100^j*sign
324+
325+ scale /= 2
326+ return(c)
327+
328+}
329+
330+# workhorse for real & imag
331+define realimag_(c,f) {
332+ auto x, sign, cs, ci, j, os
333+ os = scale
334+ cs = scale(c);ci = length(c)-cs-1
335+
336+ sign = read_digit(c,ci)
337+
338+ scale = cs/2
339+
340+ x = 0
341+ for(j=-scale;j*2<ci;j++) {
342+ x += 10^j * read_digit(c,j*2+f)
343+ }
344+
345+ if ( (sign == 2+f) || (sign == 4) )x*=-1
346+
347+ scale = os
348+ return(x)
349+}
350+
351+# Extract the real part of a generated complex number
352+define real(c) {
353+ return realimag_(c,1)
354+}
355+
356+# Extract the imaginary part of a generated complex number
357+define imag(c) {
358+ return realimag_(c,0)
359+}
360+
361+# Print a generated complex number
362+# (normal print functions don't work on complex numbers)
363+define printc(c) {
364+ auto r,i
365+ r = real(c)/1
366+ i = imag(c)/1
367+ print r
368+ if(i<0){print " -";i*=-1}else{print " +"}
369+ print " i*"
370+ return(i)
371+}
372+
373+## Basic math - ordinary bc operators don't work on these special
374+## interleaved complex numbers
375+
376+# Add two complex numbers and return another complex number
377+define cadd(c1,c2) {
378+ return makecomplex(real(c1)+real(c2),imag(c1)+imag(c2))
379+}
380+
381+# Subtract a complex number from another, returning complex
382+define csub(c1,c2) {
383+ return makecomplex(real(c1)-real(c2),imag(c1)-imag(c2))
384+}
385+
386+# Multiply two complex, return complex
387+define cmul(c1,c2) {
388+ auto r1,r2,i1,i2;
389+ r1 = real(c1)
390+ i1 = imag(c1)
391+ r2 = real(c2)
392+ i2 = imag(c2)
393+ return makecomplex(r1*r2-i1*i2, r2*i1+r1*i2)
394+}
395+
396+# Calculates the absolute value of a complex number
397+# NB: returns a standard bc number and not a new fangled 'complex'
398+define cabs(c) {
399+ return sqrt(real(c)^2+imag(c)^2)
400+}
401+
402+# Returns the complex conjugate of a complex number
403+# ...is faster than the hashed out return line
404+define cconj(c) {
405+ #return makecomplex(real(c), -imag(c))
406+ return c - 2*10^(length(c)-scale(c)-1)
407+}
408+
409+# Divide one complex by another, returning a third
410+define cdiv(c1,c2) {
411+ auto r1,r2,i1,i2,aa;
412+ r1 = real(c1)
413+ r2 = real(c2)
414+ i1 = imag(c1)
415+ i2 = imag(c2)
416+ aa = r2^2+i2^2
417+ return makecomplex((r1*r2+i1*i2)/aa, (r2*i1-r1*i2)/aa)
418+}
419+
420+# Negate a complex number (i.e. multiply by -1)
421+# ...faster than cmul(c, makecomplex(-1,0))
422+# ...or the hashed out return line
423+define cneg(c) {
424+ auto sign, t;
425+ # return makecomplex(-real(c), -imag(c))
426+ t = length(c) - scale(c) - 1
427+ sign = read_digit(c,t)
428+ c += (5-2*sign)*(10^t)
429+ return(c)
430+}
431+
432+## Basic Math II - fast[er] squaring, square roots and integer power
433+
434+# Square a complex number
435+define csquare(c) {
436+ auto r,i
437+ r = real(c)
438+ i = imag(c)
439+ return makecomplex(r^2-i^2, 2*r*i)
440+}
441+
442+# Find the primary square root of a complex number
443+define csqrt(c) {
444+ auto r,i, sr, si
445+ r = real(c)
446+ i = imag(c)
447+ if(i==0){if(r>=0){
448+ return(makecomplex(sqrt(r),0))
449+ } else {
450+ return(makecomplex(0,sqrt(-r)))
451+ }}
452+ sr = sqrt((sqrt(r^2+i^2)+r)/2)
453+ si = i/sr/2
454+ return makecomplex(sr,si)
455+}
456+
457+# subfunctions for use by cintpow
458+define mod(n,m) {auto s;s=scale;scale=0;n%=m;scale=s;return(n)}
459+define int(n) {auto s;s=scale;scale=0;n/=1;scale=s;return(n)}
460+
461+# Raise a complex number [c] to an integer power [n]
462+# NB: n is a standard bc number and not a complex
463+define cintpow(c, n) {
464+ auto x
465+ x = makecomplex(1,0)
466+ n = int(n)
467+ if(n<0)return( cdiv(x,cintpow(c,-n)) )
468+ if(n==0)return( x )
469+ if(n==1)return( c )
470+ if(n==2)return( csquare(c) )
471+ while(n) {
472+ if(mod(n,2)) {
473+ x = cmul(x, c) ; n -= 1
474+ #print "DEBUG: cintpow - odd_ (",n,")\n"
475+ } else {
476+ c = csquare(c) ; n = int(n/2)
477+ #print "DEBUG: cintpow - even (",n,")\n"
478+ }
479+ }
480+ return x
481+}
482
483=== added file 'bc-lib/electric.bc'
484--- bc-lib/electric.bc 1970-01-01 00:00:00 +0000
485+++ bc-lib/electric.bc 2009-03-21 02:39:06 +0000
486@@ -0,0 +1,314 @@
487+### electric.bc - a large number of functions for use with GNU BC
488+
489+## Not to be regarded as suitable for any purpose
490+## Not guaranteed to return correct answers
491+#
492+# jmj(at)bowlfr.org
493+
494+# Calcule du nombre pi
495+
496+# (loi d'OHM) : Calcul de la résistance d'un circuit connaissant la tension appliquée et l'intensité du courant.
497+# Énoncé : La résistance, exprimée en ohms, s'obtient en divisant la tension exprimée en volts par l'intensité du courant, exprimée en ampères.
498+# La formule peut aussi s'appliquer en exprimant la tension et l'intensité du courant respectivement en sous-multiples (ou en multiples) du volt et de l'ampère ; en particulier, on peut exprimer la tension en millivolts (mV ; 1 mV = 0,001 V) et l'intensité du courant en milliampères (mA ; 1 mA = 0,001 A) et ainsi la résistance sera exprimée en ohms.
499+# Cependant, si la tension est exprimée en volts et l'intensité du courant en milliampères (mA ; 1 mA = 0,001 A), la résistance sera exprimée en kilo-ohms (kW ; 1 kW = 1 000 W). On exprime quelquefois la tension en volts et l'intensité du courant en microampères (µA ; 1 µA = 0,000 001 A) ; dans ce cas, en appliquant la formule, on obtient la résistance exprimée en mégohms (MW ; 1 MW = 1 000 000 W).
500+# r = résistance en W (ohm)
501+# v = tension en V (volt)
502+# i = intensité du courant en A (ampère)
503+
504+define r2vi(v,i) {
505+ return( v/i )
506+}
507+
508+# (loi d'OHM) : Calcul de la tension appliquée à un circuit connaissant la résistance et l'intensité du courant
509+# La formule peut aussi s'appliquer en exprimant la résistance et l'intensité du courant respectivement avec les sous-multiples de l'ohm et de l'ampère ; on peut en particulier exprimer la résistance en kilo-ohms (kW) et l'intensité du courant en milliampères (mA) et ainsi la tension sera exprimée en volts.
510+# v = tension en V (volt)
511+# r = résistance en W (ohm)
512+# i = intensité du courant en A (ampère)
513+
514+define v2ri(r,i) {
515+ return( r*i )
516+}
517+
518+# (loi d'OHM) : Calcul de l'intensité du courant connaissant la résistance et la tension appliquée au circuit.
519+# La formule peut aussi s'appliquer en exprimant la tension et la résistance respectivement avec les sous-multiples et multiples du volt et de l'ohm. Par exemple, la tension peut être exprimée en millivolts (mV) et la résistance en ohm (W) : dans ce cas, l'intensité du courant sera exprimée en milliampères (mA) ; de même, la tension peut être exprimée en volts et la résistance en kilo-ohms (kW), l'intensité du courant sera exprimée en milliampères.
520+# i = intensité du courant en A (ampère)
521+# v = tension en V (volt)
522+# r = résistance en W (ohm)
523+
524+define i2vr(v,r) {
525+ return( v/r )
526+}
527+
528+# Calcul de la résistance totale d'un circuit formé de plusieurs résistances reliées en série connaissant leurs valeurs.
529+# Énoncé : La valeur de la résistance équivalente à plusieurs résistances reliées en série est obtenue en faisant la somme des valeurs de chaque résistance.
530+# Les valeurs de résistance doivent toutes être exprimées dans la même unité de mesure (W, kW, MW).
531+# rt = résistance équivalente totale
532+# r1 = valeur de la première résistance
533+# r2 = valeur de la deuxième résistance
534+# rn = valeur de la dernière résistance
535+
536+define r2rt () {
537+ i = 1
538+ print " Veuillez entrer les valeurs des résistances.\n"
539+ print " Sortez en mettant 0 comme valeur.\n\n"
540+ while (1) {
541+ "Valeur R" ; i
542+ r = read()
543+ if (r == 0) break;
544+ rt += r
545+ i += 1
546+ }
547+ return( rt )
548+}
549+
550+# Calcul de la conductance totale d'un circuit formé de plusieurs résistances reliées en parallèle connaissant leurs valeurs.
551+# Énoncé : La conductance résultante de l'association de plusieurs résistances reliées en parallèle s'obtient en faisant la somme des conductances de chaque résistance, (voir théorie 2 dans la rubrique sommaire électronique).
552+# Les valeurs des conductances doivent toutes être exprimées dans la même unité de mesure (S ou mS ou µS).
553+# gt = conductance totale
554+# g1 = conductance de la première résistance
555+# g2 = conductance de la deuxième résistance
556+# gn = conductance de la dernière résistance
557+
558+define c2ct () {
559+ i = 1
560+ print " Veuillez entrer les valeurs des conductances.\n"
561+ print " Sortez en mettant 0 comme valeur.\n\n"
562+ while (1) {
563+ "Valeur G" ; i
564+ g = read()
565+ if (g == 0) break;
566+ gt += g
567+ i += 1
568+ }
569+ return( gt )
570+}
571+
572+# Calcul de la résistance équivalente à plusieurs résistances reliées en parallèle connaissant leurs valeurs.
573+# Énoncé : La résistance équivalente à plusieurs résistances reliées en parallèle s'obtient en effectuant les calculs en trois temps : d'abord on calcule la conductance de chaque résistance ; puis on calcule la conductance totale des résistances en parallèle ; enfin, on calcule la résistance équivalente, c'est-à-dire la résistance qui correspond à la conductance totale
574+# Les valeurs de résistance doivent toutes être exprimées dans la même unité de mesure (W, kW, MW).
575+# re = résistance équivalente
576+# r1 = valeur de la première résistance
577+# r2 = valeur de la deuxième résistance
578+# rn = valeur de la dernière résistance
579+
580+define r2re () {
581+ i = 1
582+ print " Veuillez entrer les valeurs des résistances.\n"
583+ print " Sortez en mettant 0 comme valeur.\n\n"
584+ while (1) {
585+ "Valeur R" ; i
586+ r = read()
587+ if (r == 0) break;
588+ rtemp += 1/r
589+ i += 1
590+ rt = 1/rtemp
591+ }
592+ return( rt )
593+}
594+
595+# Calcul de la résistance équivalente de deux résistances reliées en parallèle connaissant leurs valeurs.
596+# Énoncé : La résistance équivalente de deux résistances reliées en parallèle s'obtient en multipliant les valeurs des deux résistances et en divisant le tout par la somme de ces deux valeurs.
597+# Il est possible d'utiliser r2re () pour ce calcul
598+# Les valeurs de résistance doivent toutes être exprimées dans la même unité de mesure (W, kW, MW).
599+# re = résistance équivalente
600+# r1 = valeur de la première résistance
601+# r2 = valeur de la deuxième résistance
602+
603+define re (r1,r2) {
604+ return( (r1*r2)/( r1+r2) )
605+}
606+
607+# Calcul de la valeur de la résistance à mettre en parallèle avec une autre résistance de valeur connue pour obtenir une résistance équivalente donnée.
608+# Énoncé : La valeur de la résistance à mettre en parallèle avec une autre résistance de valeur connue pour obtenir une résistance équivalente donnée se calcule en multipliant la valeur de la résistance connue par la résistance équivalente, le tout divisé par la différence de ces deux valeurs.
609+# ri = résistance inconnue
610+# r = valeur de la résistance disponible
611+# re = résistance équivalente que l'on veut obtenir
612+# Les valeurs de résistance doivent toutes être exprimées dans la même unité de mesure (W, kW, MW).
613+
614+define ri (r,re) {
615+ return( (r*re)/(r-re) )
616+}
617+
618+# Calcul de la résistance équivalente à deux ou plusieurs résistances de même valeur reliées en parallèle.
619+# Énoncé : La résistance équivalente de deux ou plusieurs résistances de valeurs égales reliées en parallèle s'obtient en divisant la valeur par le nombre de résistances.
620+# La résistance équivalente sera exprimée dans la même unité de mesure que celle utilisée pour exprimer la valeur des résistances.
621+# Re = résistance équivalente
622+# R = valeur des résistances
623+# n = nombre de résistances
624+
625+define rn2re (r,n) {
626+ return( (r/n )
627+}
628+
629+# Calcul de la force électromotrice (f.e.m.) obtenue en reliant en série deux ou plusieurs piles, connaissant la force électromotrice de chaque pile.
630+# Énoncé : En mettant en série deux ou plusieurs piles, on obtient une force électromotrice égale à la somme des forces électromotrices de chaque pile.
631+# Les forces électromotrices doivent toutes être exprimées dans la même unité de mesure.
632+# et = force électromotrice totale
633+# e1 = force électromotrice de la première pile
634+# e2 = force électromotrice de la deuxième pile
635+# en = force électromotrice de la dernière pile.
636+
637+define en2et () {
638+ i = 1
639+ print " Veuillez entrer les valeurs des forces électromotrices.\n"
640+ print " Sortez en mettant 0 comme valeur.\n\n"
641+ while (1) {
642+ "Valeur E" ; i
643+ e = read()
644+ if (r == 0) break;
645+ et += e
646+ i += 1
647+ }
648+ return( et )
649+}
650+
651+# Calcul de la résistance interne d'une pile connaissant sa f.e.m. (tension à vide lorsqu'elle ne fournit aucun courant) et la tension en charge lorsqu'elle fournit un courant donné.
652+# Énoncé : La résistance interne d'une pile est donnée par la différence entre la f.e.m. et la tension en charge, le tout divisé par le courant fourni.
653+# ri = résistance interne en W (ohm)
654+# e = f.e.m. "tension à vide" en V (volt)
655+# v = tension en charge en V (volt)
656+# i = courant fournit en A (ampère).
657+
658+define evi2ri (e,v,i) {
659+ return( (e-v)/i )
660+}
661+
662+# Calcul de la puissance électrique d'un appareil connaissant la tension appliquée et l'intensité du courant absorbé.
663+# Énoncé : La puissance électrique, exprimée en watts, s'obtient en multipliant la tension, exprimée en volts, par l'intensité du courant exprimée en ampères.
664+# p = puissance électrique en W (watt)
665+# v = tension en V (volt)
666+# i = intensité du courant en A (ampère)
667+
668+define vi2p (v,i) {
669+ return( v*i )
670+}
671+
672+# Calcul de l'intensité du courant absorbé par un appareil connaissant la tension appliquée et sa puissance électrique.
673+# i = intensité du courant en A (ampère)
674+# p = puissance électrique en W (watt)
675+# v = tension appliquée en V (volt)
676+
677+define pv2i (p,v) {
678+ return( p/v )
679+}
680+
681+# Calcul de la tension appliquée à un appareil connaissant l'intensité du courant absorbé et la puissance électrique de celui-ci.
682+# v = tension appliquée en V (volt)
683+# p = puissance électrique en W (watt)
684+# i = intensité du courant absorbé en A (ampère).
685+
686+define pi2v (p,i) {
687+ return( p/i )
688+}
689+
690+# Calcul de l'énergie électrique consommée par un appareil connaissant sa puissance électrique et sa durée de fonctionnement.
691+# Énoncé : L'énergie consommée par un appareil, exprimée en watts-secondes, s'obtient en multipliant la puissance de l'appareil, exprimée en watts, par le temps de fonctionnement exprimé en secondes.
692+# Le watt-seconde (W.s), unité de mesure utilisée pour exprimer la quantité d'énergie électrique consommée, est équivalent à 1 joule (J), unité de mesure de l'énergie et du travail mécanique.
693+# 1 W.s = 1 J
694+# En pratique, pour indiquer la consommation domestique et industrielle de l'énergie électrique, on utilise un multiple du watt-seconde, c'est-à-dire le kilowatt-heure (kW.h).
695+# 1 kW.h = 3 600 000 W.s ; 1 W.s = 1 / 3 600 000 kW.h
696+# Le résultat de l'exemple précédent peut s'exprimer en kW.h au moyen de l'équivalence :
697+# 180 000 W.s = 180 000 / 3 600 000 kW.h = 0,05 kW.h
698+# Si dans la formule 88, la puissance est exprimée en kilowatts (kW ; 1 kW = 1000 W) et le temps en heures (h, 1 h = 60 mn = 3 600 s), l'énergie consommée sera exprimée en kilowatts-heures.
699+# Par exemple, si P = 800 W = 0,8 kW et "t" = 30 mn = 0,5 h, l'énergie consommée sera de :
700+# W = 0,8 x 0,5 = 0,4 kW.h
701+# w = énergie consommée en W.s (watt-seconde)
702+# p = puissance électrique en W (watt)
703+# t = temps en s (secondes)
704+
705+define pt2w (p,t) {
706+ return( p*t )
707+}
708+
709+# Calcul de la puissance électrique dissipée par effet Joule dans une résistance connaissant l'intensité du courant et la valeur de cette résistance.
710+# Énoncé : La puissance électrique, exprimée en watts, dissipée dans une résistance s'obtient en multipliant la résistance, exprimée en ohms, par le carré du courant qui la traverse, exprimé en ampères.
711+# p = puissance électrique en W (watt)
712+# r = résistance en W (ohm)
713+# i = intensité du courant en A (ampère)
714+
715+define ri2p (r,i) {
716+ return( r*i^2 )
717+}
718+
719+# Calcul de la valeur d'une résistance connaissant la puissance électrique dissipée et l'intensité du courant qui la traverse.
720+# r = résistance en W (ohm)
721+# p = puissance dissipée en W (watt)
722+# i = intensité du courant en A (ampère)
723+
724+define pi2r (p,i) {
725+ return( p/i^2 )
726+}
727+
728+# Calcul de l'intensité du courant qui parcourt une résistance connaissant la puissance électrique dissipée et la valeur de cette résistance.
729+# Énoncé : L'intensité du courant, exprimée en ampères, qui parcourt une résistance, s'obtient en divisant la puissance dissipée exprimée en watts par la résistance exprimée en ohms, et en extrayant la racine carrée du quotient obtenu.
730+# i = intensité du courant en A (ampère)
731+# p = puissance dissipée en W (watt)
732+# r = résistance en W (ohm)
733+
734+define pr2i (p,r) {
735+ return( sqrt(p/r) )
736+}
737+
738+# Calcul de la puissance électrique dissipée par effet Joule dans une résistance connaissant la tension appliquée et la valeur de cette résistance.
739+# Énoncé : La puissance électrique, exprimée en watts, dissipée dans une résistance s'obtient en divisant le carré de la tension, exprimée en volts, par la résistance exprimée en ohms.
740+# p = puissance électrique en W (watt)
741+# v = tension en V (volt)
742+# r = résistance en W (ohm)
743+
744+define vr2p (v,r) {
745+ return( v^2/r )
746+}
747+
748+# Calcul de la valeur d'une résistance connaissant la puissance électrique dissipée et la tension appliquée.
749+# r = résistance en W (ohm)
750+# v = tension appliquée en V (volt)
751+# p = puissance électrique dissipée en W (watt)
752+
753+define vp2r (v,p) {
754+ return( v^2/p )
755+}
756+
757+# Calcul de la tension appliquée à une résistance connaissant la puissance électrique dissipée et la valeur de cette résistance.
758+# Énoncé : La tension appliquée à une résistance, exprimée en volts, s'obtient en multipliant la puissance dissipée exprimée en watts, par la résistance exprimée en ohms et en extrayant la racine carrée du produit obtenu.
759+# v = tension en V (volt)
760+# p = puissance dissipée en W (watt)
761+# r = résistance en W (ohm)
762+
763+define pr2v (p,r) {
764+ return( sqrt(p*r) )
765+}
766+
767+# Calcul de la quantité de chaleur obtenue en transformant par effet Joule une quantité d'énergie électrique donnée (pour le calcul de l'énergie électrique, voir la formule pt2w(p,t).
768+# Énoncé : La quantité de chaleur exprimée en kilo-calories produite dans une résistance par effet Joule s'obtient en multipliant l'énergie électrique exprimée en watts-secondes (Joule) dissipée dans cette résistance par le nombre 0,000238.
769+# qc = 0,000238 W (Valeur approchée par défaut)
770+# qc = quantité de chaleur en kilo-calories
771+# w = énergie électrique en W.s (watt-seconde) ou bien en J (Joule)
772+# Données : Énergie électrique dissipée par une résistance W = 0,5 kW.h (kilowatt-heure) = 3 600 000 x 0,5 = 1 800 000 W.s (pour l'équivalence entre le kilowatt-heure et le watt-seconde, voir l'observation qui suit la formule pt2w(p,t).
773+# Quantité de chaleur produite par la résistance : Qc = 0,000 238 x 1 800 000 = 428,4 kcal.
774+
775+define qcw2qc (w) {
776+ # convertion en W.s ou Joule
777+ return( qc )
778+}
779+
780+# Calcul de la résistance à chaud d'un conducteur connaissant l'augmentation de température du matériau et la résistance du conducteur à la température ambiante (20° C).
781+# Énoncé : En augmentant la température d'un conducteur, on augmente sa résistance électrique. Pour le calcul de la résistance à chaud d'un conducteur, il faut compléter l'énoncé de la façon suivante : la résistance à chaud, exprimée en ohms, s'obtient en additionnant la résistance du conducteur à la température ambiante (20° C) avec le produit du coefficient de température du matériau, de la résistance à la température ambiante et de l'augmentation de température exprimée en degrés Celsius.
782+# Les coefficients de température des principaux matériaux conducteurs sont reportés dans la dernière colonne de droite du tableau III (nous reportons la même figure 1 ci-dessous).
783+# rt = résistance à chaud (à la température t) en W (ohm)
784+# r20 = résistance à froid (à la température de 20° C) en W (ohm)
785+# a = coefficient de température du matériau
786+# t = température du conducteur chaud en °C (degrés Celsius)
787+# t - 20 = augmentation de température en °C (degrés Celsius)
788+# Données relatives à un conducteur de tungstène : R20 = 30 W (résistance à froid du conducteur) ; = 0,0045 (coefficient de température du tungstène) ; t = 320° C (température du conducteur).
789+# Augmentation de température du conducteur : t - 20 = 320 - 20 = 300° C
790+# Résistance à chaud du conducteur : Rt = 30 + 0,0045 x 30 x 300 = 30 + 40,5 = 70,5 W
791+
792+
793+# Calcul de la résistance par mètre à froid (à 20° C) d'un conducteur connaissant sa section et la résistivité du matériau.
794+# (R / m) = (p / S)
795+# R / m = résistance par mètre en W / m (ohm par mètre)
796+# p = résistivité en µW.m (microhm-mètre)
797+# S = section en mm2
798+# * Cette formule est tirée de la formule 64 ("voir formulaire mathématiques 2 - 1ère partie") en donnant à la longueur du conducteur la valeur de 1 mètre *.
799+
800+
801
802=== added file 'bc-lib/funcs.bc'
803--- bc-lib/funcs.bc 1970-01-01 00:00:00 +0000
804+++ bc-lib/funcs.bc 2009-03-21 02:39:06 +0000
805@@ -0,0 +1,447 @@
806+### Funcs.BC - a large number of functions for use with GNU BC
807+
808+## Not to be regarded as suitable for any purpose
809+## Not guaranteed to return correct answers
810+
811+scale=100;
812+define pi(){return(a(1)*4)} ; pi = pi()
813+ e = e(1)
814+define phi(){return((1+sqrt(5))/2)} ; phi = phi()
815+define psi(){return((1-sqrt(5))/2)} ; psi = psi()
816+
817+### Reset base to ten
818+define resetbase() { ibase=1+1+1+1+1+1+1+1+1+1;obase=ibase; }
819+d0=0;d1=1;d2=2;d3=3;d4=4;d5=5;d6=6;d7=7;d8=8;d9=9
820+d10=10;d11=11;d12=12;d13=13;d14=14;d15=15;d16=16;d17=17;d18=18;d19=19
821+d20=20;d21=21;d22=22;d23=23;d24=24;d25=25;d26=26;d27=27;d28=28;d29=29
822+d30=30;d31=31;d32=32;d33=33;d34=34;d35=35;d36=36;d37=37;d38=38;d39=39
823+
824+### Integer and Rounding
825+
826+# Round to next integer nearest 0: -1.99 -> 1, 0.99 -> 0
827+define int(x) { auto os;os=scale;scale=0;x/=1;scale=os;return(x) }
828+
829+# Round down to integer below x
830+define floor(x) { auto xx;xx=int(x);if(xx>x)xx-=1;return(xx) }
831+
832+# Round up to integer above x
833+define ceil(x) { return(-floor(-x)) }
834+
835+# Fractional part of x: 12.345 -> 0.345
836+define frac(x) { return(x-floor(x)) }
837+
838+# Absolute value of x
839+define abs(x) { if(x<0)return(-x)else return(x) }
840+
841+# Sign of x
842+define sgn(x) { if(x<0)return(-1)else if(x>0)return(1);return(0) }
843+
844+# Round x up to next multiple of y
845+define rup(x,y) { return(y*ceil(x/y)) }
846+
847+# Round x down to previous multiple of y
848+define rdn(x,y) { return(y*floor(x/y)) }
849+
850+# Round x to the nearest multiple of y
851+define rnd(x,y) { return(y*floor(x/y+.5)) }
852+
853+# Find the remainder of x/y
854+define rem(x,y) { return(x-rdn(x,y)) }
855+
856+# Greatest common divisor of x and y
857+define gcd(x,y) { auto r;while(y>0){r=rem(x,y);x=y;y=r};return(x) }
858+
859+# Lowest common multiple of x and y
860+define lcm(x,y) { return (x*y/gcd(x,y)) }
861+
862+# Output function - prints a and b as a fraction in smallest terms
863+define sft(a,b) { #smallest fractional terms
864+ auto c,d,e
865+ c=gcd(a,b);
866+ d=int(a/c);
867+ e=int(b/c);
868+ print a,"/",b," = ",d,"/",e,"\n";
869+ return(d/e)
870+}
871+
872+### Exponential / Logs / Trig
873+
874+# Raise x to the y-th power
875+define pow(x,y) {
876+ auto yy;yy=int(y)
877+ if (y==yy) { return(x^yy); }
878+ return( e(y*l(x)) );
879+}
880+
881+# y-th root of x [ x^(1/y) ]
882+define root(x,y) {
883+ auto iy,iyy;
884+ iy=1/y;iyy=int(iy)
885+ if (iy==iyy) { return(x^iyy); }
886+ scale+=5;y=e(l(x)/y);scale-=5;y=y/1+10^-scale
887+ x=int(y);if(x==y)y=x
888+ return( y );
889+}
890+
891+# workhorse function for powrem
892+define powrem_(x,y,m) {
893+ auto r, y2;
894+ if(y==0)return(1)
895+ if(y==1){if(x<m){return(x)}else{return(x-m*(x/m))}}
896+ y2=y/2
897+ r=powrem_(x,y2,m); if(r>=m)r%=m
898+ r^=2 ; if(r>=m)r%=m
899+ if(2*y2!=y){r*=x ; if(r>=m)r%=m}
900+ return( r )
901+}
902+
903+# Raise x to the y-th power, modulo m
904+define powrem(x,y,m) {
905+ auto os,r;
906+ os=scale;scale=0
907+ if(x<0){
908+ print "powrem: base is negative, rectifying...\n"
909+ x*=-1
910+ }
911+ if(x!=x/1){
912+ print "powrem: base is not an integer, truncating...\n"
913+ x/=1
914+ }
915+ if(y<0){
916+ print "powrem: exponent is negative, rectifying...\n"
917+ y*=-1
918+ }
919+ if(y!=y/1){
920+ print "powrem: exponent is not an integer, truncating...\n"
921+ y/=1
922+ }
923+ if(m<0){
924+ print "powrem: modulus is negative, rectifying...\n"
925+ m*=-1
926+ }
927+ if(m!=m/1){
928+ print "powrem: modulus is not an integer, truncating...\n"
929+ m/=1
930+ }
931+ if(m==0){
932+ print "powrem: modulus is zero - doing full calculation!\n"
933+ return x^y
934+ }
935+ r=powrem_(x,y,m)
936+ scale=os
937+ return( r )
938+}
939+
940+# Logarithm of y, base x: log(2, 32) = 5 because 2^5 = 32
941+define log(x,y) { return(l(y)/l(x)) }
942+
943+# Sine
944+#efine s(x) { predefined }
945+# Cosine
946+#efine c(x) { predefined }
947+# Tangent
948+define t(x) { auto c;c=c(x);if(c==0)c+=10^-scale;return(s(x)/c) }
949+
950+# Secant
951+define sc(x) { return(1/s(x)) }
952+# Cosecant
953+define cs(x) { return(1/c(x)) }
954+# Cotangent
955+define ct(x) { auto s;s=s(x);if(s==0)s+=10^-scale;return(c(x)/s) }
956+
957+# Arcsine
958+define as(x) { if(abs(x)==1)return(2*a(1)*x)else return( a(x/sqrt(1-x^2)) ) }
959+# Arccosine
960+define ac(x) { return 2*a(1)-as(x) }
961+
962+# Arctangent (one argument)
963+#efine a(x) { single argument arctangent is predefined }
964+
965+# Arctangent (two arguments)
966+define at(x,y) {
967+ auto p;
968+ if(x==0&&y==0)return(0)
969+ p=(1-sgn(y))*2*a(1)*(2*(x>=0)-1)
970+ if(x==0||y==0)return(p)
971+ return(p+a(x/y))
972+}
973+
974+# Arcsecant
975+define asc(x) { return( a(x/sqrt(x^2-1)) ) }
976+# Arccosecant
977+define acs(x) { return( asc(x)+2*a(1)*(sgn(x)-1) ) }
978+# Arctangent
979+define act(x) { return( a(x)+2*a(1) ) }
980+
981+# Hyperbolic Sine
982+define sh(x) { auto t;t=e(x);return((t-1/t)/2) }
983+# Hyperbolic Cosine
984+define ch(x) { auto t;t=e(x);return((t+1/t)/2) }
985+# Hyperbolic Tangent
986+define th(x) { auto t;t=e(2*x)-1;return(t/(t+2)) }
987+
988+# Hyperbolic Secant
989+define sch(x) { return(1/ch(x)) }
990+# Hyperbolic Cosecant
991+define csh(x) { return(1/sh(x)) }
992+# Hyperbolic Cotangent
993+define cth(x) { return(1/th(x)) }
994+
995+# Hyperbolic Arcsine
996+define ash(x) { return( l(x+sqrt(x^2+1)) ) }
997+# Hyperbolic Arccosine
998+define ach(x) { return( l(x+sqrt(x^2-1)) ) }
999+# Hyperbolic Arctangent
1000+define ath(x) { return( l((1+x)/(1-x))/2 ) }
1001+
1002+# Hyperbolic Arcsecant
1003+define asch(x) { return( l((sqrt(1-x^2)+1)/x) ) }
1004+# Hyperbolic Arccosecant
1005+define acsh(x) { return( l((sqrt(1+x^2)*sgn(x)+1)/x) ) }
1006+# Hyperbolic Arccotangent
1007+define acth(x) { return( l((x+1)/(x-1))/2 ) }
1008+
1009+define grad(x) { auto s;s=s(a(1)-x);if(s==0)s+=10^-scale;return(s(x)/s) }
1010+define agrad(x) { return( x/sqrt(1+x^2) ) }
1011+define chord(x) { return( 2* s(x/2) ) }
1012+define achord(x) { return( 2*as(x/2) ) }
1013+
1014+# Length of the diagonal vector (0,0)-(x,y) [pythagoras]
1015+define pyth(x,y) { return(sqrt(x^2+y^2)) }
1016+define pyth3(x,y,z) { return(sqrt(x^2+y^2+z^2)) }
1017+
1018+### Triangular numbers
1019+
1020+# xth triangular number
1021+define tri(x) {
1022+ auto xx
1023+ x=x*(x+1)/2;xx=int(x)
1024+ if(x==xx)return(xx)
1025+ return(x)
1026+}
1027+
1028+# 'triangular root' of x
1029+define trirt(x) {
1030+ auto xx
1031+ x=(sqrt(1+8*x)-1)/2;xx=int(x)
1032+ if(x==xx)return(xx)
1033+ return(x)
1034+}
1035+
1036+# Workhorse for following 2 functions
1037+define tri_step_(t,s) {
1038+ auto tt
1039+ t=t+(1+s*sqrt(1+8*t))/2;tt=int(t)
1040+ if(tt==t)return(tt)
1041+ return(t)
1042+}
1043+
1044+# Turn tri(x) into tri(x+1) without knowing x
1045+define tri_succ(t) {
1046+ return(tri_step_(t,0+1))
1047+}
1048+
1049+# Turn tri(x) into tri(x-1) without knowing x
1050+define tri_pred(t) {
1051+ return(tri_step_(t,0-1))
1052+}
1053+
1054+### Fibonacci
1055+
1056+# n-th Fibonacci number
1057+define fib(x){
1058+ auto a,b,c,os
1059+ os=scale;scale=0;x/=1
1060+ a=0;b=1;c=1
1061+ if(x<0){scale=os;return(fib(-x)*((-1)^(1-x)))}
1062+ if(x==0){scale=os;return(0)}
1063+ while(--x){
1064+ c=a+b;a=b;b=c
1065+ }
1066+ scale=os;return(c)
1067+}
1068+
1069+define fibf(n) { return( (pow(phi,n)-pow(psi,n))/sqrt(5) ) }
1070+
1071+# n-th Lucas number
1072+define luc(x){
1073+ auto a,b,c,os
1074+ os=scale;scale=0;x/=1
1075+ a=2;b=1;c=3
1076+ if(x<0){scale=os;return(luc(-x)*((-1)^(-x)))}
1077+ if(x==0){scale=os;return(2)}
1078+ if(x==1){scale=os;return(1)}
1079+ while(--x){
1080+ c=a+b;a=b;b=c
1081+ }
1082+ scale=os;return(c)
1083+}
1084+
1085+define lucf(n) { return( pow(phi,n)+pow(psi,n) ) }
1086+
1087+### Factorials
1088+
1089+# x!
1090+define factorial(x) {
1091+ auto i,xx
1092+ if(x<0)return(0)
1093+ if(x<2)return(1)
1094+ xx=1;for(i=x;i>=1;i--)xx*=i
1095+ return(xx)
1096+}
1097+
1098+# Gosper's approximation to the natural log of n!
1099+define gosper(n) { return( n*(l(n)-1) + ( l(2*n+1/3)+ l(pi()) )/2 ) }
1100+
1101+# Gosper's approximation to n!
1102+define gfactorial(n) { return ceil(e(gosper(n))) }
1103+
1104+# logarithm of x!
1105+define lnfactorial(x) {
1106+ auto i,xx,max
1107+ if(x<0)return(-10^100)
1108+ if(x<2)return(0)
1109+ max=2500 # Arbitrary large value
1110+ if(x<max)return( l(factorial(x)) )
1111+ xx=l(factorial(max))
1112+ for(i=x;i>max;i--)xx+=l(i)
1113+ return(xx)
1114+}
1115+
1116+# Number of permutations of r items from a group of n
1117+define permutation(n,r) {
1118+ auto i,p
1119+ if(n<0||r<0||r>n)return(0)
1120+ p=1;for(i=n;i>n-r;i--)p*=i
1121+ return(p)
1122+}
1123+
1124+# Number of combinations of r items from a group of n
1125+define combination(n,r) {
1126+ if(n<0||r<0||r>n)return(0)
1127+ if(2*r>n)r=n-r
1128+ return( permutation(n,r)/factorial(r) )
1129+}
1130+
1131+# y-th factorial of x: x!_y
1132+define multifactorial(y,x) {
1133+ auto i,xx;
1134+ xx=1;for(i=x;i>=1;i-=y)xx*=i
1135+ return(xx);
1136+}
1137+
1138+### Digit related functions
1139+
1140+# The base ten number created by appending the base ten numbers
1141+# from 1 to x, i.e. 1, 12, 123, ..., 123456789101112, etc.
1142+define dc(x) {
1143+ if (x<=0) return(0);
1144+ return(x+dc(x-1)*10^int(1+log(10,x)));
1145+}
1146+
1147+# Sum of digits in a number: 1235 -> 11. ibase modifies base used
1148+define digitsum(x) {
1149+ auto os,xx,t;
1150+ os=scale;scale=0
1151+ t=0;while(x>=1){ xx=x/ibase;t+=x-xx*ibase;x=xx }
1152+ scale=os
1153+ return(t)
1154+}
1155+
1156+# Product of digits+1 less #digits
1157+# e.g. 235 -> (2+1)(3+1)(5+1)-3 = 3*4*6 - 3 = 69
1158+define digitprod(x) {
1159+ auto os,xx,t,c;
1160+ os=scale;scale=0
1161+ t=1;c=0;while(x>=1){ xx=x/ibase;t*=(x-xx*ibase+1);x=xx;c+=1 }
1162+ scale=os
1163+ return(t-c)
1164+}
1165+
1166+# Number of digits in the base 'ibase' representation of x
1167+define digits(x) {
1168+ auto os,c;
1169+ if(x<0)return(digits(-x))
1170+ if(x==0||x==1)return(1)
1171+ os=scale;scale=10
1172+ c=ceil(log(ibase,x))-3;if(c<0)c=0
1173+ scale=0
1174+ x/=ibase^c
1175+ while(x){c+=1;x/=ibase}
1176+ scale=os
1177+ return(c)
1178+}
1179+
1180+# Reverse the digits in x (use ibase)
1181+define reverse(x) {
1182+ auto os,y;
1183+ os=scale;scale=0
1184+ y = 0
1185+ while(x) {
1186+ y = ibase*y + x%ibase
1187+ x /= ibase
1188+ }
1189+ scale=os
1190+ return(y)
1191+}
1192+
1193+### Formatted output
1194+
1195+# workhorse function for expnot and engnot
1196+define ezznot_(x,f,z) {
1197+ auto m, e, ms, es, l10, os, ns;
1198+ os=scale ; ns=scale+6 ; scale=ns
1199+ ms = 1
1200+ if (x < 0) { ms = -1 ; x *= -1 }
1201+
1202+ l10 = l(10)
1203+ x = l(x)/l10
1204+ es = 1
1205+ if (x < 0) { es = -1 ; x *= -1 }
1206+
1207+ scale=0
1208+ e = x/z; if(es==-1)e+=1 ; e*=z
1209+ scale=ns
1210+ m = e(l10*es*(x-e)) + 10^-(os+1)
1211+ scale=f
1212+ m /= 1
1213+ scale=os
1214+
1215+ print ms*m,"*10^"
1216+ return(es*e)
1217+}
1218+
1219+# Exponential notation - display x in the form a.bbbbbbb*10^cc
1220+define expnot(x,f) { # f = sig. fig.
1221+ return ezznot_(x,f,1)
1222+}
1223+
1224+# Engineering notation - display x in the form a.bbbbbbb*10^cc
1225+# where cc is a multiple of 3
1226+define engnot(x,f) { # f = sig. fig.
1227+ return ezznot_(x,f,3)
1228+}
1229+
1230+# Truncate trailing zeroes from a scaled number
1231+define trunc(x) {
1232+ auto os,i;os=scale
1233+ for(i=0;i<=os;i++){
1234+ scale=i
1235+ if(x==x/1){
1236+ x=x/1
1237+ scale=os
1238+ return(x)
1239+ }
1240+ }
1241+}
1242+
1243+# Print an integer in a field width
1244+define intprint(n, w){ # w is field width
1245+ auto os,m,i;
1246+ os=scale;scale=0;n/=1
1247+ m=n;w+=(m!=0);if(m<0){m*=-1;w-=1}
1248+ for(;m>0;w--){m/=10}
1249+ for(i=1;i<w;i++)print " "
1250+ scale=os;return(n)
1251+}
1252+#print "funcs.bc loaded\n";
1253
1254=== added file 'bc-lib/geometrie.bc'
1255--- bc-lib/geometrie.bc 1970-01-01 00:00:00 +0000
1256+++ bc-lib/geometrie.bc 2009-03-21 02:39:06 +0000
1257@@ -0,0 +1,262 @@
1258+### geometrie.bc - a large number of functions for use with GNU BC
1259+
1260+## Not to be regarded as suitable for any purpose
1261+## Not guaranteed to return correct answers
1262+#
1263+# jmj(at)bowlfr.org
1264+
1265+# Calcule du nombre pi
1266+define pi() {
1267+ scale=20
1268+ return(4*a(1))
1269+}
1270+
1271+# Calcul de la surface d'un triangle connaissant les valeurs de la base et de la hauteur
1272+# b = base
1273+# h = hauteur
1274+define s2tq(b,h) {
1275+ return( (b*h)/2 )
1276+}
1277+
1278+# Calcul de la surface d'un triangle équilatéral, "triangle ayant trois côtés égaux" connaissant la longueur du côté
1279+# c = côté
1280+define s2te(c) {
1281+ return( (sqrt(3)/4)*c^2 )
1282+}
1283+
1284+# Calcul de la surface d'un triangle isocèle "triangle ayant deux côtés égaux" connaissant la valeur des côtés égaux et de la base.
1285+# b = base
1286+# c = côté
1287+define s2ti(b,c) {
1288+ return( b/2*sqrt(c^2-(b/2)^2))
1289+}
1290+
1291+# Calcul de la surface d'un triangle scalène "triangle ayant trois côtés inégaux" connaissant la longueur des côtés.
1292+# a = côté
1293+# b = côté
1294+# c = côté
1295+define s2ts(a,b,c) {
1296+ p = (a+b+c)/2
1297+ return( sqrt(p*(p-a)*(p-b)*(p-c)))
1298+}
1299+
1300+# Calcul de l'hypoténuse d'un triangle rectangle connaissant les deux autres côtés (le triangle rectangle est un triangle ayant un angle de 90° ; l'hypoténuse est le plus grand côté, les deux autres côtés forment l'angle de 90°)
1301+# a = côté
1302+# b = côté
1303+define h2tr(a,b) {
1304+ return( sqrt(a^2+b^2))
1305+}
1306+
1307+# Calcul d'un côté d'un triangle rectangle connaissant les longueurs de l'hypoténuse et de l'autre côté
1308+# h = hypoténuse
1309+# b = côté
1310+define c2tr(h,b) {
1311+ return( sqrt(h^2-b^2))
1312+}
1313+
1314+# Calcul de la surface d'un triangle rectangle connaissant les deux côtés de l'angle droit.
1315+# a = côté
1316+# b = côté
1317+define s2tr(a,b) {
1318+ return( (a*b)/2 )
1319+}
1320+
1321+# Calcul de la diagonale d'un carré connaissant la longueur du côté
1322+# c = côté
1323+define d2ca(c) {
1324+ return( sqrt(2)*c )
1325+}
1326+
1327+# Calcul de la surface d'un carré connaissant la longueur du côté
1328+# c = côté
1329+define s2cac(c) {
1330+ return( c^2 )
1331+}
1332+
1333+# Calcul de la surface d'un carré connaissant la longueur de la diagonale
1334+# d = diagonale
1335+define s2cad(d) {
1336+ return( d^2/2 )
1337+}
1338+
1339+# Calcul de la diagonale d'un rectangle connaissant les valeurs de la base et de la hauteur
1340+# b = base
1341+# h = hauteur
1342+define d2r(b,h) {
1343+ return( sqrt(b^2+h^2) )
1344+}
1345+
1346+# Calcul de la surface d'un rectangle ou parallélogramme connaissant les valeurs de la base et de la hauteur.
1347+# b = base
1348+# h = hauteur
1349+define s2r(b,h) {
1350+ return( b*h )
1351+}
1352+
1353+# Calcul de la surface d'un losange connaissant la longueur des diagonales
1354+# (le losange est un quadrilatère ayant quatre côtés égaux et des angles adjacents inégaux).
1355+# d1 = grande diagonale
1356+# d2 = petite diagonale
1357+
1358+define s2l(d1,d2) {
1359+ return( (d1*d2)/2 )
1360+}
1361+
1362+# Calcul de la surface d'un trapèze connaissant les valeurs des deux bases et de la hauteur.
1363+# h = hauteur
1364+# b1 = petite base
1365+# b2 = grande base
1366+
1367+define s2t(h,b1,b2) {
1368+ return( (h*(b1+b2))/2 )
1369+}
1370+
1371+# Calcul de la surface d'un pentagone régulier connaissant la longueur des côtés (le pentagone régulier est un polygone ayant cinq côtés égaux et cinq angles égaux)
1372+# c = côté
1373+define s2p(c) {
1374+ return( 1.72*c^2 )
1375+}
1376+
1377+# Calcul de la surface d'un hexagone régulier connaissant la longueur d'un côté (l'hexagone régulier est un polygone ayant six côtés égaux et six angles internes égaux)
1378+# c = côté
1379+define s2h(c) {
1380+ return( 2.6*c^2 )
1381+}
1382+
1383+# Calcul du périmètre d'un cercle (circonférence) connaissant la valeur du diamètre
1384+# d = diamètre
1385+define p2c(d) {
1386+ return( pi()*d )
1387+}
1388+
1389+# Calcul de la surface d'un cercle connaissant la valeur du diamètre
1390+# d = diamètre
1391+define s2c(d) {
1392+ return( pi()/4*d^2 )
1393+}
1394+
1395+# Calcul de la longueur d'un arc de cercle connaissant la valeur de l'angle au centre et la longueur du rayon
1396+# v = valeur de l'angle au centre (en degré)
1397+# r = rayon
1398+define l2ac(v,r) {
1399+ return( pi()/180*v*r )
1400+}
1401+
1402+# Calcul de la surface d'un secteur circulaire connaissant la valeur de l'angle au centre et la longueur du rayon (un secteur circulaire est la surface plane délimitée par un arc de cercle et deux rayons).
1403+# v = valeur de l'angle au centre (en degré)grand diamètre
1404+# r = rayon
1405+define s2sc(v,r) {
1406+ return( pi()/360*v*r^2 )
1407+}
1408+
1409+# Calcul de la surface d'une couronne circulaire connaissant la valeur des deux diamètres (une couronne circulaire est la surface plane comprise entre deux circonférences concentriques).
1410+# d1 = grand diamètre
1411+# d2 = petit diamètre
1412+define s2cc(d1,d2) {
1413+ return( pi()/4*(d1^2 - d2^2) )
1414+}
1415+
1416+# Calcul de la surface d'un segment de parabole connaissant la valeur de la base et de la hauteur (on appelle segment de parabole la surface plane comprise entre un arc de parabole et la corde sous-tendue entre les extrémités de l'arc).
1417+# b = base
1418+# h = hauteur
1419+define s2sp(b,h) {
1420+ return( 2/3*b*h )
1421+}
1422+
1423+# Calcul de la surface d'une ellipse connaissant la longueur des deux axes
1424+# a = grand axe
1425+# b = petit axe
1426+define s2e(a,b) {
1427+ return( pi()/4*a*b )
1428+}
1429+
1430+# Calcul de la longueur d'une hélice connaissant le nombre de spires, les valeurs du diamètre et de la hauteur
1431+# n = nombre de spires
1432+# d = diamètre
1433+# h = hauteur
1434+define l2h(n,d,h) {
1435+ return( sqrt(pi()^2*n^2*d^2+h^2) )
1436+}
1437+
1438+# Calcul du volume d'un cube connaissant la longueur de l'arête
1439+# a = arête
1440+define v2cu(a) {
1441+ return( a^3 )
1442+}
1443+
1444+# Calcul d'une diagonale d'un cube
1445+# d = diagonale
1446+# a = arête
1447+define d2cu(a) {
1448+ return( a*sqrt(3) )
1449+}
1450+
1451+# Calcul du volume d'un parallélépipède connaissant les valeurs de la longueur et de la largeur de la base, et la hauteur
1452+# a = longueur de la base
1453+# b = largeur de la base
1454+# h = hauteur
1455+define v2pa(a,b,h) {
1456+ return( a*b*h )
1457+}
1458+
1459+# Calcul du volume d'un cylindre connaissant les valeurs du diamètre et de la hauteur
1460+# d = diamètre
1461+# h = hauteur
1462+define v2cy(d,h) {
1463+ return( pi()/4*d^2*h )
1464+}
1465+
1466+# Calcul du volume d'un cylindre creux connaissant les valeurs des deux diamètres et de la hauteur
1467+# d1 = diamètre externe
1468+# d2 = diamètre interne
1469+# h = hauteur
1470+define v2cyc(d1,d2,h) {
1471+ return( pi()/4*(d1^2-d2^2)*h )
1472+}
1473+
1474+# Calcul du volume d'un anneau à section carrée connaissant les valeurs des diamètres externes et internes
1475+# d1 = diamètre externe
1476+# d2 = diamètre interne
1477+define v2a(d1,d2) {
1478+ return( pi()/8*(d1-d2)^2*(d1+d2) )
1479+}
1480+
1481+# Calcul du volume d'un tore (anneau à section circulaire) connaissant la valeur du diamètre extérieur et celle du diamètre de la section de l'anneau
1482+# d1 = diamètre externe
1483+# d2 = diamètre interne
1484+define v2t(d1,d2) {
1485+ return( pi()^2/4*(d1-d2)*d2^2 )
1486+}
1487+
1488+# Calcul de la surface d'une sphère connaissant la valeur du diamètre
1489+# d = diamètre
1490+define s2s(d) {
1491+ return( pi()*d^2 )
1492+}
1493+
1494+# Calcul du volume d'une sphère connaissant la valeur du diamètre
1495+# d = diamètre
1496+define v2s(d) {
1497+ return( pi()/6*d^3 )
1498+}
1499+
1500+# Calcul de la surface d'une calotte sphérique connaissant les valeurs du diamètre du contour et de la hauteur
1501+# d = diamètre du contour
1502+# h = hauteur de la calotte
1503+define s2cs(d,h) {
1504+ return( pi()/4*(d^2+4*h^2 ) )
1505+}
1506+
1507+# Calcul du volume d'une calotte sphérique connaissant la valeur du diamètre de la base et de la hauteur
1508+# d = diamètre de la base
1509+# h = hauteur de la calotte
1510+define v2cs(d,h) {
1511+ return( pi()*h^2*((3*d^2+4*h^2)/(24*h)) )
1512+}
1513+
1514+# Calcul du volume d'une paraboloïde connaissant la valeur du diamètre de la base et de la hauteur
1515+# d = diamètre de la base
1516+# h = hauteur de la paraboloïde
1517+define v2pb(d,h) {
1518+ return( pi()/8*d^2*h )
1519+}
1520
1521=== added file 'bc-lib/intdiff.bc'
1522--- bc-lib/intdiff.bc 1970-01-01 00:00:00 +0000
1523+++ bc-lib/intdiff.bc 2009-03-21 02:39:06 +0000
1524@@ -0,0 +1,78 @@
1525+### IntDiff.BC - numeric differentiation and integration of
1526+### a single variable function with GNU bc
1527+
1528+## Not guaranteed for any particular purpose
1529+
1530+define f(x) { return x^2 } # example - redefine in later code/bc session
1531+
1532+# Numerically differentiate the global function f(x)
1533+define dfxdx(x) {
1534+ auto eps, os;
1535+ os = scale
1536+ scale = 0
1537+ eps = os/2
1538+ scale = os
1539+ eps = 10^-eps
1540+ return (f(x+eps)-f(x-eps))/(2*eps)
1541+}
1542+
1543+# New global variable like 'scale' - determines accuracy of numerical
1544+# integration at the expense of time. Don't set above 15 unless this
1545+# is running on a really fast machine!
1546+depth = 10
1547+
1548+# Numerically integrate the global function f(x) between x = a and x = b
1549+define ifxdx(a,b) {
1550+ auto h,s,t
1551+ h = 2^-depth
1552+ s = (b - a) * h
1553+ for(i=a+s;i<b;i+=s)t+=f(i);t+=(f(a)+f(b))/2;t*=s
1554+ return t
1555+}
1556+
1557+# glai - guess limit at infinity
1558+# Assumes p, q and r are 3 consecutive convergents to a limit and
1559+# attempts to extrapolate precisely what that limit is after an infinite
1560+# number of iterations.
1561+
1562+# 0 = glai returns function result only
1563+# 1 = glai commentates on interesting convergents
1564+glaitalk = 1
1565+
1566+define glai(p,q,r) {
1567+ auto m,n
1568+ m = q^2-p*r
1569+ n = 2*q-p-r
1570+ if(n==0)if(m==0){
1571+ if(glaitalk)print "glai: Constant series detected\n"
1572+ return p
1573+ }else{
1574+ if(glaitalk)print "glai: Arithmetic progression detected: limit is infinite\n"
1575+ return 1/0
1576+ }
1577+ if(m==0){
1578+ if(glaitalk)print "glai: Geometric progression detected: limit wraps to zero!\n"
1579+ return 0
1580+ }
1581+ return m/n
1582+}
1583+
1584+# Examples:
1585+# glai(x,x+1,x+2) causes a division by zero error as the limit of
1586+# an arithmetic progression is infinite
1587+# glai(a*k,a*k^2,a*k^3) returns zero! The limit of a geometric
1588+# progression in p-adics is precisely that,
1589+# and somehow this simple function 'knows'.
1590+# glai(63.9, 63.99, 63.999) returns 64 - correctly predicting the
1591+# limit of the sequence.
1592+
1593+# Run consecutive convergents to the ifxdx function through glai()
1594+# attaining "better" accuracy with slightly fewer calculations
1595+define ifxdx_g(a,b) {
1596+ auto p,q,r
1597+ depth -= 3 ; p = ifxdx(a,b)
1598+ depth += 1 ; q = ifxdx(a,b)
1599+ depth += 1 ; r = ifxdx(a,b)
1600+ depth += 1
1601+ return glai(p,q,r)
1602+}
1603
1604=== added file 'bc-lib/logic.bc'
1605--- bc-lib/logic.bc 1970-01-01 00:00:00 +0000
1606+++ bc-lib/logic.bc 2009-03-21 02:39:06 +0000
1607@@ -0,0 +1,103 @@
1608+### Logic.BC - Do bitwise functions with GNU bc
1609+
1610+# Perform a bitwise logical AND of x and y
1611+define and(x,y) {
1612+ auto z, t, xx, yy, os;
1613+ os=scale;scale=0
1614+ z=0;t=1;for(;x||y;){
1615+ xx=x/2;yy=y/2
1616+ z+=t*((x-2*xx)&&(y-2*yy))
1617+ t*=2;x=xx;y=yy
1618+ }
1619+ scale=os;return (z)
1620+}
1621+
1622+# Perform a bitwise logical OR of x and y
1623+define or(x,y) {
1624+ auto z, t, xx, yy, os;
1625+ os=scale;scale=0
1626+ z=0;t=1;for(;x||y;){
1627+ xx=x/2;yy=y/2
1628+ z+=t*((x-2*xx)||(y-2*yy))
1629+ t*=2;x=xx;y=yy
1630+ }
1631+ scale=os;return (z)
1632+}
1633+
1634+# Perform a bitwise logical EXCLUSIVE-OR of x and y
1635+define xor(x,y) {
1636+ auto z, t, xx, yy, os;
1637+ os=scale;scale=0
1638+ z=0;t=1;for(;x||y;){
1639+ xx=x/2;yy=y/2
1640+ z+=t*((x-2*xx)!=(y-2*yy))
1641+ t*=2;x=xx;y=yy
1642+ }
1643+ scale=os;return (z)
1644+}
1645+
1646+# Perform a bitwise logical NOT of x
1647+define not(x,w) {
1648+ # w is bit width
1649+ # w is implicitly expanded to the bit-width of x if
1650+ # it is not wide enough
1651+ auto z, t, xx, os;
1652+ os=scale;scale=0
1653+ z=0;t=1;for(;x||(w>0);w--){
1654+ xx=x/2;z+=t*(1-x+2*xx)
1655+ t*=2;x=xx
1656+ }
1657+ scale=os;return (z)
1658+}
1659+
1660+# Perform a bitwise logical LEFT-SHIFT of x by n places
1661+define shl(x,n,w) { # w is bit width
1662+ auto os;os=scale;scale=0
1663+ x=(x*2^n)%(2^w)
1664+ scale=os;return(x)
1665+}
1666+
1667+# Perform a bitwise logical RIGHT-SHIFT of x by n places
1668+define shr(x,n,w) { # w is bit width
1669+ auto os;os=scale;scale=0
1670+ x=(x/2^n)%(2^w)
1671+ scale=os;return(x)
1672+}
1673+
1674+# Reverse bits in x
1675+define bitrev(x,w) { # w is bit width, w = 0 -> use number if bits in x
1676+ auto os,z;os=scale;scale=0
1677+ if(w){x%=(2^w)}else{for(;2^w<x;w++){}}
1678+ z=0;for(;w;w--){z=z*2+(x%2);x/=2}
1679+ scale=os;return(z)
1680+}
1681+
1682+# Perform XOR on binary floating point representations of x and y
1683+define xorf(x,y){
1684+ auto os,t
1685+ os=scale;scale=0
1686+ t=2^(4*os)
1687+ x*=t;y*=t
1688+ scale=os
1689+ return( xor(x,y)/t )
1690+}
1691+
1692+# Perform OR on binary floating point representations of x and y
1693+define orf(x,y){
1694+ auto os,t
1695+ os=scale;scale=0
1696+ t=2^(4*os)
1697+ x*=t;y*=t
1698+ scale=os
1699+ return( or(x,y)/t )
1700+}
1701+
1702+# Perform AND on binary floating point representations of x and y
1703+define andf(x,y){
1704+ auto os,t
1705+ os=scale;scale=0
1706+ t=2^(4*os)
1707+ x*=t;y*=t
1708+ scale=os
1709+ return( and(x,y)/t )
1710+}
1711
1712=== added file 'bc-lib/pdhi.bc'
1713--- bc-lib/pdhi.bc 1970-01-01 00:00:00 +0000
1714+++ bc-lib/pdhi.bc 2009-03-21 02:39:06 +0000
1715@@ -0,0 +1,69 @@
1716+# pdhi - Pan Digital Halving Index with GNU bc
1717+# Returns how many times x must be divided by 2 before
1718+# the result contains all digits from 0 to 9 (if ibase = 10).
1719+# e.g. 3339 -> 1669.5 -> 834.75 -> 417.375 ->
1720+# 208.6875 -> 104.34375 -> 52.171875 ->
1721+# 26.0859375 -> 13.04296875, i.e. 8 times
1722+
1723+# Uses ibase as the base for divisions (usually 10)
1724+
1725+define pdhi(x) {
1726+ auto d[],xi,xf,c,r,pdhi,lim,i;
1727+ if(x<0)x*=-1
1728+ c=1;pdhi=-1;lim=int(10/ibase+3)*scale
1729+ while(c){
1730+ pdhi+=1
1731+ xi=int(x);xf=x-xi
1732+ while(xi){
1733+ r=int(xi/ibase)
1734+ d[xi-ibase*r]=1
1735+ xi=r
1736+ }
1737+ for(i=lim ; i && xf ; i--){
1738+ #while(xf){
1739+ xf*=ibase
1740+ r=int(xf)
1741+ d[r]=1
1742+ xf-=r
1743+ }
1744+ c=ibase
1745+ for(r=0;r<ibase;r++){c-=d[r];d[r]=0}
1746+ x/=2
1747+ }
1748+ return pdhi;
1749+}
1750+
1751+# pdmi(x, m) - Pan Digital Multiplying Index
1752+# Returns how many times x must be multiplied by m before
1753+# the result contains all digits from 0 to 9 (if ibase = 10).
1754+# e.g. pdmi(3339,0.5) -> 1669.5 -> 834.75 -> 417.375 ->
1755+# 208.6875 -> 104.34375 -> 52.171875 ->
1756+# 26.0859375 -> 13.04296875, i.e. 8 times
1757+
1758+# Uses ibase as the base for divisions (usually 10)
1759+
1760+define pdmi(x,m) {
1761+ auto d[],xi,xf,c,r,pdmi,lim,i;
1762+ if(x<0)x*=-1
1763+ c=1;pdmi=-1;lim=int(10/ibase+3)*scale
1764+ while(c){
1765+ pdmi+=1
1766+ xi=int(x);xf=x-xi
1767+ while(xi){
1768+ r=int(xi/ibase)
1769+ d[xi-ibase*r]=1
1770+ xi=r
1771+ }
1772+ for(i=lim ; i && xf ; i--){
1773+ #while(xf){
1774+ xf*=ibase
1775+ r=int(xf)
1776+ d[r]=1
1777+ xf-=r
1778+ }
1779+ c=ibase
1780+ for(r=0;r<ibase;r++){c-=d[r];d[r]=0}
1781+ x*=m
1782+ }
1783+ return pdmi;
1784+}
1785
1786=== added file 'bc-lib/primes.bc'
1787--- bc-lib/primes.bc 1970-01-01 00:00:00 +0000
1788+++ bc-lib/primes.bc 2009-03-21 02:39:06 +0000
1789@@ -0,0 +1,131 @@
1790+### Primes.BC - Primes and factorisation (rudimentary) with GNU bc
1791+
1792+# funcs.bc is required to run this library
1793+
1794+# Returns 2, 3, or number of form 6n[+-]1
1795+define aq(x) {
1796+ if(x<0)return(-aq(-x))
1797+ if(x<3)return(x+1)
1798+ x-=3;x+=int(x/2)
1799+ return(x+x+5)
1800+}
1801+
1802+# Inverse of the above
1803+define iaq(x) {
1804+ if(x<0)return(-iaq(-x))
1805+ if(x<4)return(x-1)
1806+ return((rem(x+3,6)+x+x)/6+1)
1807+}
1808+
1809+# Returns 2, 3, 5 or number of form 30n[+-]{1,7,11,13}
1810+define aq30(x) {
1811+ auto os, e, r, rh, s
1812+ os=scale;scale=0;x/=1
1813+ if(x<0){x=-aq30(-x);scale=os;return(x)}
1814+ if(x<4){x=x+1+x/3 ;scale=os;return(x)}
1815+ x-=3 ; e=x/8
1816+ r=x%8 ; rh=r/4
1817+ s=1-2*rh ; r=s*r+7*rh
1818+ scale=os;return( 30*(e+rh)-s*(r^2-7*r-1) )
1819+}
1820+
1821+# Inverse of the above
1822+define iaq30(x) {
1823+ auto os, e, r
1824+ os=scale;scale=0;x/=1
1825+ if(x<0){x=-iaq30(-x);scale=os;return(x)}
1826+ if(x<7){x=x-1-x/5 ;scale=os;return(x)}
1827+ e=x/30;r=x%30
1828+ #r+=(31-r)/30#not needed with bc integer division
1829+ r=r/6+(r-2)/7
1830+ scale=os;return(8*e+r +3)
1831+}
1832+
1833+# Cyrek's Approximation to the Prime Counting Function pi(x)
1834+define aprimepi(x) {
1835+ auto a,b,lx,k
1836+ lx=l(x)/l(10)
1837+ a=pow(lx,2/(3*lx))
1838+ b=1+2/(17*pow(ch(lx-e(2)),1/32))
1839+ k=1+l(a/b)
1840+ return(x*k/l(x))
1841+}
1842+
1843+# Output the prime factors of x
1844+define fac(x) {
1845+ auto os,j[],ji,p,n,c,xx;#oldscale,jump,jump-index,prime,nth,count,xx
1846+ os=scale;scale=0;x/=1
1847+ if(x<0){print"-1*\\\n";x*=-1}
1848+ if(x==0){return 0}
1849+ j[0]=4;j[1]=2;j[2]=4;j[3]=2;j[4]=4;j[5]=6;j[6]=2;j[7]=6
1850+ for(p=2;p<7;p+=p-1){
1851+ xx=x/p
1852+ c=0;while(x==p*xx){c+=1;x=xx;xx=x/p}
1853+ if(c){print p;if(c>1){print "^",c};print "*\\\n"}
1854+ }
1855+ ji=7;p=7
1856+ while(p^2<=x){
1857+ xx=x/p
1858+ c=0;while(x==p*xx){c+=1;x=xx;xx=x/p}
1859+ if(c){print p;if(c>1){print "^",c};print "*\\\n"}
1860+ p+=j[ji=(ji+1)%8]
1861+ }
1862+ if(x>1)print x,"*\\\n"
1863+ scale=os;return(1)
1864+}
1865+
1866+# Determine whether x is prime
1867+define is_prime(x) {
1868+ auto os,j[],ji,p,n,xx;#oldscale,jump,jump-index,prime,nth,xx
1869+ os=scale;scale=0;x/=1
1870+ if(x<0){return 0}
1871+ j[0]=4;j[1]=2;j[2]=4;j[3]=2;j[4]=4;j[5]=6;j[6]=2;j[7]=6
1872+ for(p=2;p<7;p+=p-1){if(x==p){return 1};xx=x/p;if(x==p*xx){return 0}}
1873+ ji=7;p=7
1874+ while(p^2<=x){
1875+ xx=x/p;if(x==p*xx){return 0}
1876+ p+=j[ji=(ji+1)%8]
1877+ }
1878+ scale=os;return(1)
1879+}
1880+
1881+# Determine whether x is y-th power-free
1882+# e.g. has_freedom(51, 2) will return whether 51 is square-free
1883+# Special case: has_freedom(x, 1) returns whether x is prime
1884+define has_freedom(x,y) {
1885+ auto os,j[],ji,p,py,n,xx;#oldscale,jump,jump-index,prime,nth,xx
1886+ os=scale;scale=0;x/=1
1887+ if(y==1){return is_prime(x)}
1888+ if(x<0||y<1){return 0}
1889+ j[0]=4;j[1]=2;j[2]=4;j[3]=2;j[4]=4;j[5]=6;j[6]=2;j[7]=6
1890+ for(p=2;p<7;p+=p-1){if(x==p){return 1};py=p^y;xx=x/py;if(x==py*xx){return 0}}
1891+ ji=7;p=7;py=p^y
1892+ while(py<=x){
1893+ xx=x/py;if(x==py*xx){return 0}
1894+ p+=j[ji=(ji+1)%8];py=p^y
1895+ }
1896+ scale=os;return(1)
1897+}
1898+
1899+# Fill the global prime[] array with g primes
1900+define genprimes(g) {
1901+ auto os,j[],ji,i,x,xx,p,p2,pflag
1902+ #os-oldscale,j[]-jumps,ji-jumpindex,i-primeindex,x-primecandidate
1903+ #xx-x/p,p-prime,pflag-is.x.prime?
1904+ os=scale;scale=0
1905+ j[0]=4;j[1]=2;j[2]=4;j[3]=2;j[4]=4;j[5]=6;j[6]=2;j[7]=6
1906+ prime[0]=1;prime[1]=2;prime[2]=3;prime[3]=5;prime[4]=7;primetop=4
1907+ if(g<100)g=100;if(g>65535)g=65535
1908+ for(i=primetop+1;i<=g;i++)prime[i]=0
1909+ ji=0;x=11
1910+ while(primetop<g){
1911+ i=2;p=prime[i];pflag=1
1912+ while(p^2<=x){
1913+ xx=x/p;if(x==p*xx){pflag=0;break}
1914+ i+=1;p=prime[i]
1915+ }
1916+ if(pflag){primetop+=1;prime[primetop]=x;}#print primetop,": ",x,"\n"}
1917+ x+=j[ji=(ji+1)%8]
1918+ }
1919+ scale=os;return(prime[primetop])
1920+}
1921
1922=== added file 'bc-lib/ss.bc'
1923--- bc-lib/ss.bc 1970-01-01 00:00:00 +0000
1924+++ bc-lib/ss.bc 2009-03-21 02:39:06 +0000
1925@@ -0,0 +1,32 @@
1926+## ss.bc calculate the value of sum[n=1..oo] x^(2^(-n))-1 with GNU bc
1927+## i.e. sum of repeated square roots, less one.
1928+
1929+define ss(x) {
1930+ auto s,t,os;
1931+ if(x<=0){print "Negative Infinity\n";-1/0}
1932+ scale+=6
1933+ s=x ; os=s+1
1934+ t=0;while(os!=s){
1935+ os=s
1936+ s=sqrt(s)
1937+ t+=s-1
1938+ }
1939+ scale-=6;t/=1
1940+ return(t);
1941+}
1942+
1943+define ss_n(n,x) { # Generalisation of ss; ss(x) == ss_n(2,x)
1944+ auto s,t,os;
1945+ if(x<=0){print "Negative Infinity\n";-1/0}
1946+ if(n==2)return(ss(x))
1947+ if(n<=1){print "Infinity\n";1/0}
1948+ scale+=6
1949+ s=x ; os=s+1
1950+ t=0;while(os!=s){
1951+ os=s
1952+ s=e(l(s)/n)
1953+ t+=s-1
1954+ }
1955+ scale-=6;t/=1
1956+ return(t);
1957+}
1958
1959=== added file 'bc-lib/thermometer.bc'
1960--- bc-lib/thermometer.bc 1970-01-01 00:00:00 +0000
1961+++ bc-lib/thermometer.bc 2009-03-21 02:39:06 +0000
1962@@ -0,0 +1,70 @@
1963+### Temperature.BC - Conversions of temperature scales to other scales for GNU bc
1964+
1965+## Celcius to... ###########################################
1966+
1967+# ... Farenheit
1968+define c2f(c) { return (c * 1.8 + 32) }
1969+
1970+# ... Kelvin
1971+define c2k(c) { return (c + 273.15) }
1972+
1973+# ... Reamur
1974+define c2re(c) { return (c * 0.8) }
1975+
1976+# ... Rankine
1977+define c2ra(c) { return (c * 1.8 + 491.67) }
1978+
1979+## Farenheit to... ########################################
1980+
1981+# ... Celcius
1982+define f2c(f) { return ((f - 32)/1.8) }
1983+
1984+# ... Kelvin
1985+define f2k(f) { return ((f + 459.67)/1.8) }
1986+
1987+# ... Reamur
1988+define f2re(f) { return ((f - 32)/2.25) }
1989+
1990+# ... Rankine
1991+define f2ra(f) { return (f + 459.67) }
1992+
1993+## Kelvin to... ###########################################
1994+
1995+# ... Celcius
1996+define k2c(k) { return (k - 273.15) }
1997+
1998+# ... Farenheit
1999+define k2f(k) { return (k * 1.8 - 459.67) }
2000+
2001+# ... Reamur
2002+define k2re(k) { return ((k - 273.15)*0.8) }
2003+
2004+# ... Rankine
2005+define k2ra(k) { return (k * 1.8) }
2006+
2007+## Reamur to... ##########################################
2008+
2009+# ... Celcius
2010+define re2c(r) { return (r / 0.8) }
2011+
2012+# ... Farenheit
2013+define re2f(r) { return (r * 2.25 + 32) }
2014+# ... Kelvin
2015+define re2k(r) { return (r / 0.8 + 273.15) }
2016+
2017+# ... Rankine
2018+define re2ra(r) { return (r * 2.25 + 491.67) }
2019+
2020+## Rankine to... #########################################
2021+
2022+# ... Celcius
2023+define ra2c(r) { return (r / 1.8 + 273.15) }
2024+
2025+# ... Farenheit
2026+define ra2f(r) { return (r - 459.67) }
2027+
2028+# ... Kelvin
2029+define ra2k(r) { return (r / 1.8) }
2030+
2031+# ... Reamur
2032+define ra2re(r) { return ((r / 1.8 + 273.15)*0.8) }

Subscribers

People subscribed via source and target branches

to all changes: