Chinese remainder theorem: Difference between revisions
Added Algol 68
Thundergnat (talk | contribs) m (Automated syntax highlighting fixup (second round - minor fixes)) |
(Added Algol 68) |
||
| (14 intermediate revisions by 8 users not shown) | |||
Line 129:
x= 23
</pre>
=={{header|AArch64 Assembly}}==
{{works with|as|Raspberry Pi 3B version Buster 64 bits <br> or android 64 bits with application Termux }}
<syntaxhighlight lang AArch64 Assembly>
/* ARM assembly AARCH64 Raspberry PI 3B */
/* program chineserem64.s */
/************************************/
/* Constantes */
/************************************/
/* for this file see task include a file in language AArch64 assembly*/
.include "../includeConstantesARM64.inc"
/*********************************/
/* Initialized data */
/*********************************/
.data
szMessResult: .asciz "Result = "
szCarriageReturn: .asciz "\n"
.align 2
arrayN: .quad 3,5,7
arrayA: .quad 2,3,2
.equ ARRAYSIZE, (. - arrayA)/8
/*********************************/
/* UnInitialized data */
/*********************************/
.bss
sZoneConv: .skip 24
/*********************************/
/* code section */
/*********************************/
.text
.global main
main:
ldr x0,qAdrarrayN // N array address
ldr x1,qAdrarrayA // A array address
mov x2,#ARRAYSIZE // array size
bl chineseremainder
ldr x1,qAdrsZoneConv
bl conversion10 // call décimal conversion
mov x0,#3
ldr x1,qAdrszMessResult
ldr x2,qAdrsZoneConv // insert conversion in message
ldr x3,qAdrszCarriageReturn
bl displayStrings // display message
100: // standard end of the program
mov x0, #0 // return code
mov x8,EXIT
svc #0 // perform the system call
qAdrszCarriageReturn: .quad szCarriageReturn
qAdrsZoneConv: .quad sZoneConv
qAdrszMessResult: .quad szMessResult
qAdrarrayA: .quad arrayA
qAdrarrayN: .quad arrayN
/******************************************************************/
/* compute chinese remainder */
/******************************************************************/
/* x0 contains n array address */
/* x1 contains a array address */
/* x2 contains array size */
chineseremainder:
stp x1,lr,[sp,-16]! // save registers
stp x2,x3,[sp,-16]! // save registers
stp x4,x5,[sp,-16]! // save registers
stp x6,x7,[sp,-16]! // save registers
stp x8,x9,[sp,-16]! // save registers
mov x4,#1 // product
mov x5,#0 // sum
mov x6,#0 // indice
1:
ldr x3,[x0,x6,lsl #3] // load a value
mul x4,x3,x4 // compute product
add x6,x6,#1
cmp x6,x2
blt 1b
mov x6,#0
mov x7,x0 // save entry
mov x8,x1
mov x9,x2
2:
mov x0,x4 // product
ldr x1,[x7,x6,lsl #3] // value of n
sdiv x2,x0,x1
mov x0,x2 // p
bl inverseModulo
mul x0,x2,x0 // = product / n * invmod
ldr x3,[x8,x6,lsl #3] // value a
madd x5,x0,x3,x5 // sum = sum + (result1 * a)
add x6,x6,#1
cmp x6,x9
blt 2b
sdiv x1,x5,x4 // divide sum by produc
msub x0,x1,x4,x5 // compute remainder
100:
ldp x8,x9,[sp],16 // restaur registers
ldp x6,x7,[sp],16 // restaur registers
ldp x4,x5,[sp],16 // restaur registers
ldp x2,x3,[sp],16 // restaur registers
ldp x1,lr,[sp],16 // restaur registers
ret
/***************************************************/
/* Calcul modulo inverse */
/***************************************************/
/* x0 cont.quad number, x1 modulo */
/* x0 return result */
inverseModulo:
stp x1,lr,[sp,-16]! // save registers
stp x2,x3,[sp,-16]! // save registers
stp x4,x5,[sp,-16]! // save registers
stp x6,x7,[sp,-16]! // save registers
mov x7,x1 // save Modulo
mov x6,x1 // A x0=B
mov x4,#1 // X
mov x5,#0 // Y
1: //
cmp x0,#0 // B = 0
beq 2f
mov x1,x0 // T = B
mov x0,x6 // A
sdiv x2,x0,x1 // A / T
msub x0,x2,x1,x0 // B and x2=Q
mov x6,x1 // A=T
mov x1,x4 // T=X
msub x4,x2,x1,x5 // X=Y-(Q*T)
mov x5,x1 // Y=T
b 1b
2:
add x7,x7,x5 // = Y + N
cmp x5,#0 // Y > 0
bge 3f
mov x0,x7
b 100f
3:
mov x0,x5
100:
ldp x6,x7,[sp],16 // restaur registers
ldp x4,x5,[sp],16 // restaur registers
ldp x2,x3,[sp],16 // restaur registers
ldp x1,lr,[sp],16 // restaur registers
ret
/***************************************************/
/* display multi strings */
/***************************************************/
/* x0 contains number strings address */
/* x1 address string1 */
/* x2 address string2 */
/* x3 address string3 */
/* other address on the stack */
/* thinck to add number other address * 4 to add to the stack */
displayStrings: // INFO: affichageStrings
stp x1,lr,[sp,-16]! // save registers
stp x2,x3,[sp,-16]! // save registers
stp x4,x5,[sp,-16]! // save registers
add fp,sp,#48 // save paraméters address (6 registers saved * 8 bytes)
mov x4,x0 // save strings number
cmp x4,#0 // 0 string -> end
ble 100f
mov x0,x1 // string 1
bl affichageMess
cmp x4,#1 // number > 1
ble 100f
mov x0,x2
bl affichageMess
cmp x4,#2
ble 100f
mov x0,x3
bl affichageMess
cmp x4,#3
ble 100f
mov x3,#3
sub x2,x4,#4
1: // loop extract address string on stack
ldr x0,[fp,x2,lsl #3]
bl affichageMess
subs x2,x2,#1
bge 1b
100:
ldp x4,x5,[sp],16 // restaur registers
ldp x2,x3,[sp],16 // restaur registers
ldp x1,lr,[sp],16 // restaur registers
ret
/***************************************************/
/* ROUTINES INCLUDE */
/***************************************************/
/* for this file see task include a file in language AArch64 assembly */
.include "../includeARM64.inc"
</syntaxhighlight>
{{Output}}
<pre>
Result = 23
</pre>
=={{header|Action!}}==
<syntaxhighlight lang="action!">INT FUNC MulInv(INT a,b)
Line 206 ⟶ 401:
Ada.Text_IO.Put_Line(Integer'Image(Sum mod Prod));
end Chin_Rema;</syntaxhighlight>
=={{header|ALGOL 68}}==
{{Trans|C|with some error messages}}
<syntaxhighlight lang="algol68">
BEGIN # Chinese remainder theorewm - translated from the C sample #
PROC mul inv = ( INT a in, b in )INT:
IF b in = 1
THEN 1
ELSE
INT b0 = b in;
INT a := a in, b := b in, x0 := 0, x1 := 1;
WHILE a > 1 DO
IF b = 0 THEN
print( ( "Numbers not pairwise coprime", newline ) ); stop
FI;
INT q = a OVER b;
INT t;
t := b; b := a MOD b; a := t;
t := x0; x0 := x1 - q * x0; x1 := t
OD;
IF x1 < 0 THEN x1 + b0 ELSE x1 FI
FI # mul inv # ;
PROC chinese remainder = ( []INT n, a )INT:
IF LWB n /= LWB a OR UPB n /= UPB a OR ( UPB a - LWB a ) + 1 < 1
THEN print( ( "Array bounds mismatch or empty arrays", newline ) ); stop
ELSE
INT prod := 1, sum := 0;
FOR i FROM LWB n TO UPB n DO prod *:= n[ i ] OD;
IF prod = 0 THEN
print( ( "Numbers not pairwise coprime", newline ) ); stop
FI;
FOR i FROM LWB n TO UPB n DO
INT p = prod OVER n[ i ];
sum +:= a[ i ] * mul inv( p, n[ i ] ) * p
OD;
sum MOD prod
FI # chinese remainder # ;
print( ( whole( chinese remainder( ( 3, 5, 7 ), ( 2, 3, 2 ) ), 0 ), newline ) )
END
</syntaxhighlight>
{{out}}
<pre>
23
</pre>
=={{header|ARM Assembly}}==
{{works with|as|Raspberry Pi <br> or android 32 bits with application Termux}}
<syntaxhighlight lang ARM Assembly>
/* ARM assembly Raspberry PI or android with termux */
/* program chineserem.s */
/* REMARK 1 : this program use routines in a include file
see task Include a file language arm assembly
for the routine affichageMess conversion10
see at end of this program the instruction include */
/* for constantes see task include a file in arm assembly */
/************************************/
/* Constantes */
/************************************/
.include "../constantes.inc"
/*********************************/
/* Initialized data */
/*********************************/
.data
szMessResult: .asciz "Result = "
szCarriageReturn: .asciz "\n"
.align 2
arrayN: .int 3,5,7
arrayA: .int 2,3,2
.equ ARRAYSIZE, (. - arrayA)/4
/*********************************/
/* UnInitialized data */
/*********************************/
.bss
sZoneConv: .skip 24
/*********************************/
/* code section */
/*********************************/
.text
.global main
main:
ldr r0,iAdrarrayN @ N array address
ldr r1,iAdrarrayA @ A array address
mov r2,#ARRAYSIZE @ array size
bl chineseremainder
ldr r1,iAdrsZoneConv
bl conversion10 @ call décimal conversion
mov r0,#3
ldr r1,iAdrszMessResult
ldr r2,iAdrsZoneConv @ insert conversion in message
ldr r3,iAdrszCarriageReturn
bl displayStrings @ display message
100: @ standard end of the program
mov r0, #0 @ return code
mov r7, #EXIT @ request to exit program
svc #0 @ perform the system call
iAdrszCarriageReturn: .int szCarriageReturn
iAdrsZoneConv: .int sZoneConv
iAdrszMessResult: .int szMessResult
iAdrarrayA: .int arrayA
iAdrarrayN: .int arrayN
/******************************************************************/
/* compute chinese remainder */
/******************************************************************/
/* r0 contains n array address */
/* r1 contains a array address */
/* r2 contains array size */
chineseremainder:
push {r1-r9,lr} @ save registers
mov r4,#1 @ product
mov r5,#0 @ sum
mov r6,#0 @ indice
1:
ldr r3,[r0,r6,lsl #2] @ load a value
mul r4,r3,r4 @ compute product
add r6,#1
cmp r6,r2
blt 1b
mov r6,#0
mov r7,r0 @ save entry
mov r8,r1
mov r9,r2
2:
mov r0,r4 @ product
ldr r1,[r7,r6,lsl #2] @ value of n
bl division
mov r0,r2 @ p
bl inverseModulo
mul r0,r2,r0 @ = product / n * invmod
ldr r3,[r8,r6,lsl #2] @ value a
mla r5,r0,r3,r5 @ sum = sum + (result1 * a)
add r6,#1
cmp r6,r9
blt 2b
mov r0,r5 @ sum
mov r1,r4 @ product
bl division
mov r0,r3
100:
pop {r1-r9,pc} @ restaur registers
/***************************************************/
/* Calcul modulo inverse */
/***************************************************/
/* r0 containt number, r1 modulo */
/* x0 return result */
inverseModulo:
push {r1-r7,lr} @ save registers
mov r7,r1 // save Modulo
mov r6,r1 // A r0=B
mov r4,#1 // X
mov r5,#0 // Y
1: //
cmp r0,#0 // B = 0
beq 2f
mov r1,r0 // T = B
mov r0,r6 // A
bl division // A / T
mov r0,r3 // B and r2=Q
mov r6,r1 // A=T
mov r1,r4 // T=X
mls r4,r2,r1,r5 // X=Y-(Q*T)
mov r5,r1 // Y=T
b 1b
2:
add r7,r7,r5 // = Y + N
cmp r5,#0 // Y > 0
bge 3f
mov r0,r7
b 100f
3:
mov r0,r5
100:
pop {r1-r7,pc}
/***************************************************/
/* display multi strings */
/***************************************************/
/* r0 contains number strings address */
/* r1 address string1 */
/* r2 address string2 */
/* r3 address string3 */
/* other address on the stack */
/* thinck to add number other address * 4 to add to the stack */
displayStrings: @ INFO: affichageStrings
push {r1-r4,fp,lr} @ save des registres
add fp,sp,#24 @ save paraméters address (6 registers saved * 4 bytes)
mov r4,r0 @ save strings number
cmp r4,#0 @ 0 string -> end
ble 100f
mov r0,r1 @ string 1
bl affichageMess
cmp r4,#1 @ number > 1
ble 100f
mov r0,r2
bl affichageMess
cmp r4,#2
ble 100f
mov r0,r3
bl affichageMess
cmp r4,#3
ble 100f
mov r3,#3
sub r2,r4,#4
1: @ loop extract address string on stack
ldr r0,[fp,r2,lsl #2]
bl affichageMess
subs r2,#1
bge 1b
100:
pop {r1-r4,fp,pc}
/***************************************************/
/* ROUTINES INCLUDE */
/***************************************************/
.include "../affichage.inc"
</syntaxhighlight>
{{output}}
<pre>
Result = 23
</pre>
=={{header|Arturo}}==
Line 543 ⟶ 967:
23
</pre>
=={{header|
<syntaxhighlight lang="coffeescript">crt = (n,a) ->
sum = 0
Line 568 ⟶ 992:
23
</pre>
=={{header|Common Lisp}}==
Using function ''invmod'' from [[http://rosettacode.org/wiki/Modular_inverse#Common_Lisp]].
Line 768 ⟶ 1,193:
=={{header|EasyLang}}==
{{trans|C}}
<syntaxhighlight lang="text">
func mul_inv a b .
if
if x1
.
.
return x1
.
prod = 1
sum = 0
for i
prod *= n[i]
.
for i
p = prod / n[i]
.
.
n[] = [ 3 5 7 ]
a[] = [ 2 3 2 ]
print h
</syntaxhighlight>
{{out}}
<pre>23</pre>
=={{header|EchoLisp}}==
'''egcd''' - extended gcd - and '''crt-solve''' - chinese remainder theorem solve - are included in math.lib.
Line 997 ⟶ 1,425:
10 11 4 22 9 19 3 >>>crt<<< .
</pre>
=={{header|Fortran}}==
Written in Fortran-77 style (fixed form, GO TO), directly translated from the problem description.
<syntaxhighlight lang="Fortran">
* RC task: use the Chinese Remainder Theorem to solve a system of congruences.
FUNCTION crt(n, residues, moduli)
IMPLICIT INTEGER (A-Z)
DIMENSION residues(n), moduli(n)
p = product(moduli)
crt = 0
DO 10 i = 1, n
m = p/moduli(i)
CALL egcd(moduli(i), m, r, s, gcd)
IF (gcd .ne. 1) GO TO 20 ! error exit
10 crt = crt + residues(i)*s*m
crt = modulo(crt, p)
RETURN
20 crt = -1 ! will never be negative, so flag an error
END
* Compute egcd(a, b), returning x, y, g s.t.
* g = gcd(a, b) and a*x + b*y = g
*
SUBROUTINE egcd(a, b, x, y, g)
IMPLICIT INTEGER (A-Z)
g = a
u = 0
v = 1
w = b
x = 1
y = 0
1 IF (w .eq. 0) RETURN
q = g/w
u next = x - q*u
v next = y - q*v
w next = g - q*w
x = u
y = v
g = w
u = u next
v = v next
w = w next
GO TO 1
END
PROGRAM Chinese Remainder
IMPLICIT INTEGER (A-Z)
PRINT *, crt(3, [2, 3, 2], [3, 5, 7])
PRINT *, crt(3, [2, 3, 2], [3, 6, 7]) ! no solution
END
</syntaxhighlight>
{{Out}}
<pre>
23
-1
</pre>
=={{header|FreeBASIC}}==
Partial {{trans|C}}. Uses the code from [[Greatest_common_divisor#Recursive_solution]] as an include.
Line 1,034 ⟶ 1,528:
print chinese_remainder(n(), a())</syntaxhighlight>
{{out}}<pre>23</pre>
=={{header|Frink}}==
This example solves an extended version of the Chinese Remainder theorem by allowing an optional third parameter <CODE>d</CODE> which defaults to 0 and is an integer. The solution returned is the smallest solution >= d. (This optional parameter is common in many/most real-world applications of the Chinese Remainder Theorem.)
Line 1,608 ⟶ 2,103:
<syntaxhighlight lang="matlab">>> chineseRemainder([2 3 2], [3 5 7])
ans = 23</syntaxhighlight>
=={{header|Maxima}}==
<syntaxhighlight lang="maxima">
/* Function that checks pairwise coprimality */
check_pwc(lst):=block(
sublist(cartesian_product_list(makelist(i,i,length(lst)),makelist(i,i,length(lst))),lambda([x],x[1]#x[2])),
makelist([lst[%%[i][1]],lst[%%[i][2]]],i,length(%%)),
makelist(apply('gcd,%%[i]),i,length(%%)),
if length(unique(%%))=1 and first(unique(%%))=1 then true)$
/* Chinese remainder function */
c_remainder(A,N):=if check_pwc(N)=false then "chinese remainder theorem not applicable" else block(
cn:apply("*",N),
makelist(gcdex(cn/N[i],N[i]),i,1,length(N)),
makelist(A[i]*%%[i][1]*cn/N[i],i,1,length(N)),
apply("+",%%),
mod(%%,cn));
Alis:[2,3,2]$
Nlis:[3,5,7]$
c_remainder(Alis,Nlis);
</syntaxhighlight>
{{out}}
<pre>
23
</pre>
=={{header|Modula-2}}==
<syntaxhighlight lang="modula2">MODULE CRT;
Line 2,692 ⟶ 3,213:
say 'no solution found.' /*oops, announce that solution ¬ found.*/</syntaxhighlight>
{{out|output|text= is identical to the 1<sup>st</sup> REXX version.}} <br><br>
=={{header|RPL}}==
{{trans|Python}}
{{works with|Halcyon Calc|4.2.9}}
{| class="wikitable" ≪
! RPL code
! Comment
|-
|
≪
DUP ROT 1 R→C ROT 0 R→C
'''WHILE''' DUP RE '''REPEAT'''
OVER RE OVER RE / FLOOR
OVER * NEG ROT +
'''END'''
DROP C→R ROT MOD
SWAP 1 == SWAP 0 IFTE
≫ ‘<span style="color:blue">MODINV</span>’ STO
≪ → n a
≪ 0
1 1 n SIZE '''FOR''' j n j GET * '''NEXT'''
1 n SIZE '''FOR''' j
DUP n j GET / FLOOR
ROT OVER n j GET <span style="color:blue">MODINV</span> a j GET * ROT * +
SWAP '''NEXT'''
MOD
≫ ≫ ‘<span style="color:blue">NCHREM</span>’ STO
|
<span style="color:blue">MODINV</span> ''( a b → x ) with ax = 1 mod b''
3:b 2:(r,u) ← (a,1) 1:(r',u') ← (b,0)
While r' ≠ 0
q ← r // r'
(r - q*r', u - q*u')
Forget (r',u') and calculate u mod b
Test r and return zero if a and b are not co-prime
<span style="color:blue">NCHREM</span> ''( n a → remainder )''
sum = 0
prod = reduce(lambda a, b: a*b, n)
for n_i, a_i in zip(n, a):
p = prod / n_i
sum += a_i * mul_inv(p, n_i) * p
// reorder stack for next iteration
return sum % prod
|}
'''RPL 2003 version'''
{{works with|HP|49}}
≪ → a n
≪ a 1 GET n 1 GET →V2
2 a SIZE '''FOR''' j
a j GET n j GET →V2 ICHINREM
'''NEXT'''
V→ MOD
≫ ≫ '<span style="color:blue">NCHREM</span>' STO
{ 2 3 2 } { 3 5 7 } <span style="color:blue">NCHREM</span>
{{out}}
<pre>
1: 23.
</pre>
=={{header|Ruby}}==
Brute-force.
Line 2,735 ⟶ 3,320:
p chinese_remainder([10,4,9], [11,22,19]) #=> nil
</syntaxhighlight>
=={{header|Rust}}==
<syntaxhighlight lang="rust">fn egcd(a: i64, b: i64) -> (i64, i64, i64) {
Line 3,203 ⟶ 3,789:
23 = 3 (mod 5)
23 = 2 (mod 7)</pre>
=={{header|Wren}}==
{{trans|C}}
<syntaxhighlight lang="
var mulInv = Fn.new { |a, b|
if (b == 1) return 1
Line 3,242 ⟶ 3,829:
23
</pre>
=={{header|XPL0}}==
{{trans|C}}
When N are not pairwise coprime, the program aborts due to division by zero, which is one way to convey error.
<syntaxhighlight lang "XPL0">func MulInv(A, B); \Returns X where rem((A*X) / B) = 1
int A, B;
int B0, T, Q;
int X0, X1;
[B0:= B; X0:= 0; X1:= 1;
if B = 1 then return 1;
while A > 1 do
[Q:= A / B;
T:= B; B:= rem(A/B); A:= T;
T:= X0; X0:= X1 - Q*X0; X1:= T;
];
if X1 < 0 then X1:= X1 + B0;
return X1;
];
func ChineseRem(N, A, Len);
int N, A, Len;
int P, I, Prod, Sum;
[Prod:= 1; Sum:= 0;
for I:= 0 to Len-1 do Prod:= Prod*N(I);
for I:= 0 to Len-1 do
[P:= Prod / N(I);
Sum:= Sum + A(I) * MulInv(P,N(I)) * P;
];
return rem(Sum/Prod);
];
int N, A;
[N:= [ 3, 5, 7 ];
A:= [ 2, 3, 2 ];
IntOut(0, ChineseRem(N, A, 3)); CrLf(0);
]</syntaxhighlight>
{{out}}
<pre>
23
</pre>
=={{header|zkl}}==
{{trans|Go}}
| |||